Gnss atmospheric estimation with ambiguity fixing

ABSTRACT

Methods and apparatus are presented for estimating environmental parameters from GNSS signals in real time. Some embodiments estimate a float solution using a federated ionospheric filter. Some embodiments fix ambiguities for improved estimates.

CROSS REFERENCE TO RELATED APPLICATIONS

The following are related hereto and incorporated herein in theirentirety by this reference: International Patent Publication WO2010/021656 A2 dated 25 Feb. 2010 (TNL A-2526PCT); International PatentPublication WO 2010/021658 A2 dated 25 Feb. 2010 (TNL A-2525PCT);International Patent Publication WO 2010/021659 A2 dated 25 Feb. 2010(TNL A-2524PCT); International Patent Publication WO 2010/021657 A2dated 25 Feb. 2010 (TNL A-2523PCT); International Patent Publication WO2010/021660 A2 dated 25 Feb. 2010 (TNL A-2339PCT); U.S. ProvisionalApplication for Patent No. 61/189,382 filed 19 Aug. 2008 (TNL A-2339P);U.S. patent application Ser. No. 12/224,451 filed 26 Aug. 2008, UnitedStates Patent Application Publication US 2009/0027625 A1 (TNL A-1789US);International Patent Application PCT/US07/05874 filed 7 Mar. 2007,International Publication No. WO 2008/008099 A2 (TNL A-1789PCT); UnitedStates Patent Application Publication US 2009/0224969 A1 (TNL A-1743US);International Publication No. WO 2007/032947 A1 (TNL A-1743PCT); U.S.Pat. No. 7,432,853 granted 7 Oct. 2008; (TNL A-1403US); InternationalPublication Number WO 2005/045463 A1 (TNL A-1403PCT); InternationalPublication Number WO 2011/034616 A2 (TNL A-2585PCT); InternationalPatent Application Number PCT/US2011/024733 (TNL A-2633PCT); ProvisionalApplication for U.S. Pat. No. 61/337,980 filed 14 Feb. 2010 (TNLA-2633P); and Provisional Application for U.S. Pat. No. 61/396,676 filed30 May 2010 (TNL A-2751P).

TECHNICAL FIELD

The present invention relates to the field of Global NavigationSatellite Systems GNSS). More particularly, the present inventionrelates to methods and apparatus for processing of GNSS data to estimateenvironmental parameters.

BACKGROUND ART

Over the last two decades, a large number of regional GNSS referencestation networks with several thousands of receivers have beenestablished and provide network corrections for continuous real-timecentimeter level positioning (also called Real-Time Kinematicpositioning or RTK positioning). The principle of network RTKpositioning is that tropospheric and ionospheric effects are modeled inthe regional network and largely removed for the RTK positioning. ForGNSS positioning, the atmospheric effects are considered as systematicerror to be removed in the positioning process. On the other hand, thetropospheric and ionospheric delays estimated in the GNSS network arevaluable for meteorological and space weather applications.

Numerous studies have shown that the Integrated Precipitable WaterVapour (IPWV) contained in the neutral atmosphere can be retrieved usingGNSS observations with the same level of accuracy as Water VaporRadiometers(WVR). The advantage of using a GNSS network is that itprovides continuous IPWV estimation with good spatial coverage and hightemporal resolution. Nowadays, several research organizations providenear real-time and/or post-processed IPWV estimation from a GNSSreference station network with a typical delay of 30 minutes up to aday, and an accuracy of 1-2 mm. Though the accuracy is quite good, thelarge latency means the information only reflects what happened in thepast and is thus less useful for applications like Numerical WeatherPrediction (NWP).

Space weather is the study of the environmental dynamics of “geospace”:the region above the Earth's lower atmosphere including the ionosphereand the magnetosphere. Conditions on the Sun and in the solar wind,magnetosphere, and the ionosphere can affect our lives through theeffects they have on satellites, communications, navigation, and powersystems. Scientists are now studying space weather with a wide range oftools to try to learn more about the physical and chemical processestaking place in the upper atmosphere and beyond. In recent years, GNSShas become recognized as one of the premier remote sensing tools tomonitor space weather events. The signals from the GNSS satellitestravel through the ionosphere on their way to receivers on or nearEarth's surface. The free electrons populating this region of theatmosphere affect the propagation of the signals, changing their speedand direction of travel. By processing the data from a dual-frequencyGNSS receiver, it is actually possible to estimate just how manyelectrons were encountered by the signal along its travel path—the totalelectron content (TEC). TEC is the number of electrons in a column witha cross-sectional area of one square meter centered on the signal path.If a regional network of ground-based GNSS receivers is used, then a mapof TEC above the region can be constructed.

Known software programs for estimating atmospheric models from GNSSsignal data, though not in real time, include the Bernese Software fromthe Astronomical Institute, University of Bern, Version 5.0, January2007, Chapters 11 and 12, and the GAMIT/GLOBK software from theDepartment of Earth, Atmospheric and Planetary Sciences, MassachusettsInstitute of Technology, GAMIT Reference Manual Release 10.3, 1 June2009, Chapter 7.

United States Patent Publication 2009/0224969 A1 of Peter F. Kolbdescribes a Kalman-filtering method to model the ionospheric effect withcarrier phase data from a regional network in real time. As only thecarrier phase data is used, the absolute level of the TEC needs a longtime to converge.

Improved methods and apparatus for estimation of atmospheric parametersfrom GNSS data in real time are desired.

SUMMARY

Some embodiments of the invention use pseudorange and carrier phaseobservations in a federated filter approach to accelerate theconvergence of TEC for real-time estimation.

Some embodiments of the invention exploit the full accuracy of carrierphase observations using carrier phase ambiguities to achieve highaccuracy TEC estimation in real time.

Some embodiments provide a method of processing GNSS signal data toestimate environmental parameter values, comprising: obtaining GNSS datacollected at a plurality of stations distributed over a region fromsignals received from GNSS satellites over multiple epochs; obtaining asatellite differential code bias (DCB) per satellite, and estimating areceiver differential code bias per station and a value of totalelectron content (TEC) per station per satellite from the GNSS data andthe satellite differential code biases using a federated ionosphericfilter.

Some embodiments further apply a geometric filter to estimate from theGNSS data a set of values for atmospheric parameters comprising at leastone of (i) a zenith total delay per station, (ii) a zenith total delayper station and a set of tropospheric gradients per station, and (ii)aslant total delay per station per satellite; obtain meteorological datafor locations within the region; and determine values over the regionfor at least one of (1) integrated precipitable water vapor (IPWV) and(2) tropospheric slant wet delay, from the set of values for atmosphericparameters and and the meteorological data.

In some embodiments, obtaining a satellite differential code bias persatellite comprises one of (i) retrieving a satellite differential codebias per satellite from an external source, and (ii) computing asatellite differential code bias per satellite from GNSS data collectedat a network of reference stations. In some embodiments, estimating avalue of total electron content per station per satellite comprises:estimating values for ionospheric model parameters and for a stochasticionospheric delay term per satellite per station in the federatedionospheric filter, and determining the value of total electron contentper station per satellite from the estimated values for the ionosphericmodel parameters and the stochastic ionospheric delay terms.

In some embodiments, applying a federated ionospheric filter comprises:for each satellite, applying to the GNSS signal data an ionosphericsubfilter to estimate values for local states representing parametersunique to that satellite and values for common states representingparameters common to all receivers, providing values for the commonstates and associated statistical information to a master filter; andpreparing updated estimates for the local states when updated values forthe common states are provided by the master filter; and applying to thevalues for the common states and the associated statistical informationa master filter to estimate updated values for the common states, and toprovide the updated values to the ionospheric subfilters.

In some embodiments, the states unique to the satellite comprise: a setof parameters characterizing ionosphere across the region, a stochasticionospheric delay term per station per satellite, and an integerambiguity per station per satellite. In some embodiments, the statescommon to all satellites comprise a receiver differential code bias perstation. Some embodiments further comprise resetting states in themaster filter to zero with infinite variance before any observationupdate at each epoch, and then applying the decorrelated observationsfrom each subfilter to the master filter. In some embodiments, the valueof total electron content per station per satellite comprises a valuemapped to vertical. In some embodiments, obtaining meteorological datacomprises obtaining surface meteorological data for locations within theregion. In some embodiments, obtaining meteorological data comprisesobtaining radiosonde temperature data for locations within the region.In some embodiments, elapsed time between obtaining the GNSS data anddetermining values over the region for total electron content and atleast one of integrated precipitable water vapor and tropospheric slantwet delay is no more than about five seconds.

Some embodiments provide apparatus for processing GNSS signal data toestimate environmental parameter values, comprising: a federatedionospheric filter adapted to estimate a receiver differential code biasper station and a value of total electron content (TEC) per station persatellite from satellite differential code biases and the GNSS signaldata.

Some embodiments further provide: a geometric filter adapted to estimatefrom the GNSS data a set of values for atmospheric parameters comprisingat least one of : (i) a zenith total delay per station, (ii) a zenithtotal delay per station and a tropospheric gradient per station, and(iii) a slant total delay per station per satellite; an atmosphericwatch module to determine values over the region for integratedprecipitable water vapor and tropospheric slant wet delay, from theestimated values for the atmospheric parameters and from themeteorological data for the region.

Some embodiments further comprise an element to obtain a code bias persatellite by one of: retrieving a code bias per satellite from anexternal source, and computing a code bias per satellite from GNSS datacollected at a network of GNSS reference stations. Some embodimentsfurther comprise: a module to estimate values for ionospheric modelparameters and for a stochastic ionospheric delay term per satellite perstation; and a module to determine the value of total electron contentper station per satellite from the estimated values for the ionosphericmodel parameters and the stochastic ionospheric delay terms.

In some embodiments, the federated ionospheric filter comprises: anionospheric subfilter for each satellite, to estimate values for localstates representing parameters unique to that satellite and values forcommon states representing parameters common to all receivers, from theGNSS signal data and a satellite differential code bias, to providevalues for the common states and associated statistical information tomaster filter; and to prepare updated estimates for the local stateswhen updated values for the common states are provided by the masterfilter; and a master filter to estimate updated values for the commonstates from the values for the common states and the associatedstatistical information, and to provide the updated values to theionospheric subfilters.

In some embodiments, the federated ionospheric filter is operative toreset states in the master filter to zero with infinite variance beforeobservation update at each epoch, and then to apply the decorrelatedobservations from each subfilter to the master filter. In someembodiments, the meteorological watch module is adapted to use surfacemeteorological data for locations within the region as at least a subsetof the meteorological data. In some embodiments, the meteorologicalwatch module is adapted to radiosonde sounding data as at least a subsetof the meteorological data.

Some embodiments provide a method of processing GNSS signal data toestimate environmental parameter values, comprising: obtaining GNSS datacollected at a plurality of stations distributed over a region fromsignals received from GNSS satellites over multiple epochs; obtaining asatellite code bias per satellite; estimating from the GNSS data and thesatellite code biases a set of float ambiguities per station persatellite and values for atmospheric parameters comprising at least oneof: (i) a total electron content per station per satellite, (ii) azenith total delay per station, (iii) a zenith total delay per stationand a set of tropospheric gradients per station, and (iv) a slant totaldelay per station per satellite; fixing the ambiguities; and estimatingfrom the GNSS data and the fixed ambiguities corrected values for theatmospheric parameters.

Some embodiments further comprise obtaining meteorological data forlocations within the region; and determining values over the region forat least one of integrated precipitable water vapor and troposphericslant wet delay, from the corrected values for the atmosphericparameters and from the meteorological data.

In some embodiments, estimating the values for atmospheric parameterscomprises applying at least one iterative filter to the GNSS data andthe satellite code biases. In some embodiments, estimating the valuesfor atmospheric parameters comprises applying a set of factorizedfilters to the GNSS data and the satellite code biases. In someembodiments, estimating the values for atmospheric parameters comprisesapplying a federated ionospheric filter to the GNSS data and thesatellite code biases. In some embodiments, the value of total electroncontent per station per satellite comprises a value mapped to vertical.In some embodiments, obtaining meteorological data comprises obtainingsurface meteorological data for locations within the region. In someembodiments, obtaining meteorological data comprises obtainingradiosonde temperature data for locations within the region. In someembodiments, the elapsed time between obtaining the GNSS data anddetermining values over the region for the atmospheric parameters is nomore than about 5 seconds.

Some embodiments provide apparatus for processing GNSS signal data toestimate environmental parameter values, comprising: at least onerecursive filter to estimate from satellite differential code biases andfrom the GNSS signal data a set of float ambiguities per station persatellite and a set of values for atmospheric parameters comprising atleast one of: (i) a total electron content per station per satellite,(ii) a zenith total delay per station, (iii) a zenith total delay perstation and a set of tropospheric gradients per station, (iv)a slanttropospheric total delay per station per satellite, and (v) a totalelectron content per station per satellite; an ambiguity-fixing moduleto prepare fixed ambiguities from the float ambiguities; and a correctormodule to prepare from the GNSS data and the fixed ambiguities correctedvalues for the atmospheric parameters.

Some embodiments further provide a meteorological watch module todetermine values for at least one of integrated precipitable water vaporand tropospheric slant wet delay from the corrected values for theatmospheric parameters.

In some embodiments, the at least one recursive filter comprises aKalman filter. In some embodiments, the at least one recursive filtercomprises a set of factorized filters. In some embodiments, the at leastone recursive filter comprises a federated ionospheric filter. In someembodiments, the value of total electron content per station persatellite comprises a value mapped to vertical. In some embodiments, themeteorological watch module is operative to process meteorological datacomprising surface meteorological data for locations within the region.In some embodiments, the meteorological watch module is operative toprocess meteorological data comprising radiosonde sounding data forlocations within the region. In some embodiments, the elapsed timebetween obtaining the GNSS data and determining values over the regionfor integrated precipitable water vapor and ionospheric slant wet delayis no more than about 5 seconds.

Some embodiments provide a computer program product comprising: acomputer usable medium having computer readable instructions physicallyembodied therein, the computer readable instructions when executed by aprocessor enabling the processor to perform a method as describedherein. Some embodiments provide a computer program comprising a set ofinstructions which when loaded in and executed by a processor enable theprocessor to perform a method as described herein.

BRIEF DESCRIPTION OF DRAWING FIGURES

These and other aspects and features of the present invention will bemore readily understood from the embodiments described below withreference to the drawings, in which:

FIG. 1 schematically illustrates a prior-art GNSS scenario;

FIG. 2 schematically illustrates the influence of ionosphere andtropospheric on signal of a GNSS network;

FIG. 3 schematically illustrates the slant path of a GNSS signal throughthe troposphere from satellite to receiver;

FIG. 4 schematically illustrates the slant path of a GNSS signal throughthe ionosphere satellite to receiver;

FIG. 5 illustrates how ionospheric parameters describe the ionosphere ata piercepoint relative to a reference point;

FIG. 6 schematically illustrates a GNSS network useful in implementingsome embodiments of the invention;

FIG. 7 is a flow chart of a real-time process for estimating atmosphericparameter values, using a federated ionospheric filtering method inaccordance with some embodiments of the invention;

FIG. 8 schematically illustrates a system using a federated ionosphericfilter in accordance with some embodiments of the invention;

FIG. 9 is a flow chart of a federated ionospheric filter process inaccordance with some embodiments of the invention;

FIG. 10 is a schematic diagram of a federated ionospheric filter inaccordance with some embodiments of the invention;

FIG. 11 is a flow chart of a federated ionospheric filter process inaccordance with some embodiments of the invention;

FIG. 12 is a schematic diagram of a federated ionospheric filter inaccordance with some embodiments of the invention;

FIG. 13 is a flow chart of an atmospheric watch process in accordancewith some embodiments of the invention;

FIG. 14 is a schematic diagram of atmospheric watch module in accordancewith some embodiments of the invention;

FIG. 15 is a flow chart of a process 1500 for estimating environmentalparameters using fixed ambiguities in accordance with some embodimentsof the invention;

FIG. 16 is a schematic diagram of a system for estimating environmentalparameter values using fixed ambiguities in accordance with someembodiments of the invention;

FIG. 17 is a flow chart of a process 1700 for estimating environmentalparameters using fixed ambiguities in accordance with some embodimentsof the invention;

FIG. 18 is a schematic diagram of a system for estimating environmentalparameters using fixed ambiguities in accordance with some embodimentsof the invention;

FIG. 19 is a block diagram of an integrated GNSS receiver system; and

FIG. 20 is a schematic diagram of a computer system.

DETAILED DESCRIPTION Part 1: INTRODUCTION

FIG. 1 schematically illustrates a prior-art GNSS scenario 100. Receiver110 receives GNSS signals from any number m of satellites in view, suchas at 120, 130 and 140. The signals pass through the earth's ionosphere150 and through the earth's troposphere 160. Each signal has multiplecarrier frequencies, such as frequencies L1 and L2. Receiver 110determines from the signals respective pseudo-ranges, PR1, PR2, . . . ,PRM, to the satellites. Pseudo-range determinations are distorted bysignal-path variations resulting from passage of the signals through theionosphere 150 and the troposphere 160, and from multipath effects, asindicated schematically at 170. Pseudo-ranges can be determined usingthe C/A code with an error of about one meter. However, the phases ofthe L1 and L2 carriers can be measured with an accuracy of 0.01-0.05cycles (corresponding to pseudo-range errors of 2 mm to 1 cm). Phasemeasurements of the carriers are influenced by the dispersive effects ofthe ionosphere, which vary over time.

FIG. 2 schematically illustrates at 200 an ionospheric shell 205 and aportion 210 of a tropospheric shell surrounding the Earth, withground-based reference stations 215, 220, 225, . . . 230 of a networkeach receiving signals from GNSS satellites 235, 240, . . . 245. Thetroposphere has a depth of, for example zero to about 11 km.Tropospheric delay affects the signals received by each referencestation in a manner depending on atmospheric temperature, pressure andhumidity in the vicinity of the reference station, as well as theelevation of the satellite relative to the reference station. The erroris about 1 mm per meter at ground level, such that the last meter of thesignal path to the reference station gives about 1 mm of error in thetropospheric model.

Various techniques are known for modeling tropospheric path delay on thesignals. See, for example, B. HOFMANN-WELLENHOF et al., GLOBALPOSITIONING SYSTEM: THEORY AND PRACTICE, 2d Ed., 1993, section 6.3.3,pp. 98-106. Tropospheric scaling (tropo-scaling) which lumps theatmospheric parameters into one tropo-scaling parameter can beimplemented in at least three ways. A first approach is to model ZenithTotal Delay (ZTD) representing tropospheric delay in a verticaldirection relative to the reference station as a value representingrange error δr, e.g., 2.58 meters. A second approach is to model the sumof one plus a scaling factor (1+S) such that tropospheric delay in thevertical direction T′=(1+S) T, where T is a constant, e.g., 1+S=1.0238.Another approach is to model S directly, e.g., S=2.38%. In general,“tropospheric effect” is all that affects different signal frequenciesin the same way (non-dispersive).

As shown in FIG. 3, except when a satellite is directly over a referencestation, signal rays penetrate the troposphere 310 in a slant path fromsatellite to receiver such as path 320 from satellite 330 to referencestation 340. The slant path of the signal ray from a given satellite toeach reference station penetrates the troposphere at an angle α which isdifferent for each satellite in view at the station. The troposphericmapping function is thus different for eachsatellite-to-reference-station combination at each epoch. The effect ofthe different slant angles can be compensated by relating thegeometry-dependent zenith tropospheric delay Tα with ageometry-independent tropospheric delay T₉₀° (Vertical T) by a mappingfunction m(α): Tα=m(α) T₉₀°.

As shown in FIG. 4, signal rays similarly penetrate the ionosphere 410in a path from satellite to receiver such as path 420 from satellite 430to reference station 440. This slant path is explicitly accounted for bythe so-called mapping function

ƒ_(mapping)(ζ)=1/cos(ζ),   (1)

where ζ is the angle of the signal ray with the line 450 perpendicularto the ionospheric sphere through the piercepoint 460 (e.g., line 1410).The slant path of the signal ray from a given satellite to eachreference station penetrates the ionosphere at a different angle and isdifferent for each reference station. The mapping function is thusdifferent for each satellite-to-reference-station combination. Theeffect of the different slant angles can be compensated by relating thegeometry-dependent Total Electron Content (TEC) with ageometry-independent VTEC (Vertical TEC) such that

TEC/ƒ_(mapping)(ζ)=TEC cos(ζ)=VTEC.   (2)

In the example of FIG. 4, the TEC determined along slant path 460corresponds to the VTEC along the line 450 perpendicular to theionospheric sphere 410 at piercepoint 460.

With the concept of the mapping function, the ionospheric advance acrossthe network area can be written as (here the uppercase i and j are to beunderstood as exponents, not indices)

$\begin{matrix}{{I( {{\Delta \; \lambda},{\Delta\phi}} )} = {{m( {{\Delta\lambda},{\Delta\phi}} )}{( {\sum\limits_{i,{j = 0}}^{\infty}{\alpha_{i,j}{\Delta\lambda}^{i}{\Delta\phi}^{j}}} ).}}} & (3)\end{matrix}$

That is, the ionospheric advance across the network area is expressed interms of its Taylor series (or any other set of orthogonal functions,such as spherical Bessel functions). For most purposes, and asillustrated here, the expansion can be stopped at first order, and theterminology a_(1,0)=a_(λ) and a_(0,1)=a_(φ), introduced. The expressiona_(0,0)=I₀ is the ionospheric advance at the reference point, whilea_(λ) and a_(φ), are the gradients in the ionosphere in the relativecoordinates. The ionosphere at the piercepoints is therefore expressedas

I _(n) ^(m) =m _(n) ^(m)(I ₀ ^(m) +a _(λ) ^(m)Δλ_(n) ^(m) +a _(φ)^(m)Δφ_(n) ^(m)).   (4)

Thus for each satellite m in view the parameters (I₀ ^(m),a_(λ)^(m),a_(φ) ^(m)) characterize the ionosphere across the network area.Those parameters are estimated, together with the carrier-phase integerambiguity and multipath states. Generally, if the expansion of Eq. (3)is carried to k-th order, the number of states introduced for theionosphere is (k+1)(k+2)/2. The other terms of Eq.(3) (m_(n) ^(m),Δλ_(n)^(m),Δφ_(n) ^(m)) are given by the geometry of the network and theposition of satellite m.

FIG. 5 illustrates how the ionospheric parameters (I₀ ^(m),a_(λ)^(m),a_(φ) ^(m)) describe the ionosphere at a piercepoint 510 relativeto a reference point 520. The ionosphere has a TEC of I₀ ^(m) at thereference point, with a slope a_(λ) ^(m) in angular direction λ and aslope a_(φ) ^(m) in angular direction φ. In the example of FIG. 5, theTEC 530 at piercepoint 510 is the sum of a contribution 540 equal to I₀^(m), a contribution 550 based on slope a_(λ) ^(m) and the angulardistance of piercepoint 520 from reference point 520 in direction λ, anda contribution 560 based on slope an and the angular distance ofpiercepoint 510 from reference point 520 in direction φ.

While a linear treatment of the ionosphere delivers excellentavailability, an even more realistic model takes into account thethickness of the ionosphere. As is known (for example from D. BILITZA,International Reference Ionosphere 2000, RADIO SCIENCE 2 (36) 2001,261), the electron density of the ionosphere has a certain profile ƒ(h)as a function of altitude h which peaks sharply at a height between300-400 kilometers above ground. The electron content that a rayexperiences from satellite m to station n is expressed by the integral

$\begin{matrix}{{I_{n}^{m} \propto {\int_{({x^{m},y^{m},z^{m}})}^{({x_{n},y_{n},z_{n}})}{{{sf}(h)}}}},} & (5)\end{matrix}$

where s is the measure along the direct line of sight between stationand satellite. For a simple shell model, ƒ(h)=Δ(h−h₀) (Dirac Deltadistribution), and this expression reduces to the previous mappingfunction as

$\begin{matrix}{{\frac{s}{s}_{h_{0}}} = {\frac{1}{\cos \; \phi}.}} & (6)\end{matrix}$

Using suitable parameters for ƒ(h), the integral for allstation-satellite pairs can be numerically computed at each epoch. Forpractical purposes an approximation in terms of a box profile issufficient and delivers improvements over the shell model. It is assumedthat the gradients in the ionosphere do not depend on altitude. Thisassumption can easily be relaxed by adding further gradient states fordifferent altitudes. That the finite thickness of the ionosphere is animportant feature of the model can be understood by picturing the entryand exit point of the ray of a low elevation satellite, e.g., as shownin FIG. 8 of United States Patent Application Publication US2009/0224969 A1. If the thickness of the ionospheric shell is 200kilometers, the entry point and exit point might be separated by some1000 kilometers. With typical gradients of a_(λ),a_(φ)˜10⁻³ m/km, thecontributions to the calculation of ionospheric advance differ greatlyfrom entry point to exit point.

FIG. 6 schematically illustrates a system 600 useful in implementingsome embodiments of the invention. Reference stations of an optionalglobal (worldwide) tracking network, such as reference stations 605,610, . . . 615, are distributed about the Earth with the aim of havingsubstantially continuous observability of most or all GNSS satellites.The position of each reference station is known very precisely, e.g.,within less than 2 cm. Each reference station is equipped with anantenna and tracks the GNSS signals transmitted by the satellites inview at that station, such as GNSS satellites 620, 625, . . . 630. TheGNSS signals have codes modulated on each of two or more carrierfrequencies. Each reference station of the global network acquires GNSSdata representing, for each satellite in view at each epoch,carrier-phase (carrier) observations of at least two carriers, andpseudorange (code) observations of the respective codes modulated on atleast two carriers. The reference stations also obtain the broadcastnavigation message with almanac and ephemerides of the satellites fromthe satellite signals. The almanac contains the rough position of allsatellites of the GNSS, while the so-called broadcast ephemeridesprovide more precise predictions (ca. 1 m) of the satellites' positionsand the satellites' clock error (ca. 1.5 m) over specific timeintervals.

GNSS data collected at the reference stations of the optional globalnetwork are transmitted via communications channels 635 to a globalnetwork processor 640. Global network processor 640 uses the GNSS datafrom the reference stations of the global network with other informationto estimate parameters such as satellite differential code biases (DCBs)used in regional network processing. Such a global system is describedin U.S. Provisional Application for Patent No. 61/277,184 filed 19 Sep.2009 (TNL A-2585P).

Still referring to FIG. 6, reference stations of a regional (local)tracking network, such as reference stations 645, 650, . . . 655, aredistributed over a region of the Earth so as to observe GNSS satellitesvisible over the region. The position of each reference station is knownprecisely, e.g., within less than 2 cm. Each reference station isequipped with an antenna and tracks the GNSS signals transmitted by thesatellites in view at that station, such as GNS satellites 620, 625, . .. 630. Each reference station of the regional network acquires GNSS datarepresenting, for each satellite in view at each epoch, carrier-phase(carrier) observations of at least two carriers, and pseudorange (code)observations of the respective codes modulated on at least two carriers.The regional reference stations typically also obtain the broadcastnavigation message with almanac and ephemerides of the satellites fromthe satellite signals.

GNSS data collected at the reference stations of the regional networkare transmitted via communications channels 660 to a regional networkprocessor 665. Regional network processor 665 uses the GNSS data fromthe reference stations of the regional network with other information toestimate values of atmospheric parameters as described below. Thesevalues are made available for display and/or for use in one or moreweather forecasting processors 670. The data are transmitted for exampleas shown in FIG. 6 via communications channel/s 675, an uplink 680, acommunications satellite 685, and a downlink 690. Any other suitabletransmission medium may be used including but not limited to internet,radio broadcast or mobile telephone link.

PART 2: FEDERATED IONOSPHERIC FILTER

Some embodiments in accordance with the invention estimate environmentalparameters using a federated ionospheric filter.

GPS L₁ and L₂ carrier phase and pseudorange observations can beexpressed as:

$\begin{matrix}{L_{1} = {{\lambda_{1}\phi_{1}} = {\rho + T + I + {c \cdot ( {t_{r} - t^{s}} )} + b_{1}^{r} - b_{1}^{s} + {\lambda_{1}N_{1}} + v_{1}}}} & (7) \\{L_{2} = {{\lambda_{2}\phi_{2}} = {\rho + T + {\frac{\lambda_{2}^{2}}{\lambda_{1}^{2}}I} + {c \cdot ( {t_{r} - t^{s}} )} + b_{2}^{r} - b_{2}^{s} + {\lambda_{2}N_{2}} + v_{2}}}} & (8) \\{P_{1} = {\rho + T - I + {c \cdot ( {t_{r} - t^{s}} )} + B_{1}^{r} - B_{1}^{s} + ɛ_{1}}} & (9) \\{P_{2} = {\rho + T - {\frac{\lambda_{2}^{2}}{\lambda_{1}^{2}}I} + {c \cdot ( {t_{r} - t^{s}} )} + B_{2}^{r} - B_{2}^{s} + ɛ_{2}}} & (10)\end{matrix}$

where:

L₁ and L₂ are the L₁ and L₂ carrier phase observations in meters,

φ₁ and φ₂ are the L₁ and L₂ carrier phase observations in cycles,

P₁ and P₂ are the L1 and L2 pseudorange observations in meters,

ρ is the geometric range between satellite and receiver antenna phasecenter,

T is the tropospheric delay in meters,

I is the L₁ ionospheric delay in meters,

t^(s) and t_(r) are the satellite and receiver clock error,

b₁ ^(s) and b₂ ^(s) are the satellite L₁ and L₂ phase biases,

b₁ ^(r) and b₂ ^(r) are the receiver L₁ and L₂ phase biases,

B₁ ^(s) and B₂ ^(s) are the satellite L₁ and L₂ code biases,

B₁ ^(r) and B₂ ^(r) are the receiver L₁ and L₂ code biases,

N₁ and N₂ are L₁ and L₂ integer ambiguities,

ν₁ and ν₂ are L₁ and L₂ phase noise plus multipath, and

ε₁ and ε₂ are L₁ and L₂ code noise plus multipath.

Ionospheric phase and pseudorange observation mapped to L1 frequency canbe written as:

$\begin{matrix}{L_{I} = {I + {\frac{\lambda_{1}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}( {b_{2}^{r} - b_{1}^{r}} )} - {\frac{\lambda_{1}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}( {b_{2}^{s} - b_{1}^{s}} )} + N_{I} + v_{I}}} & (11) \\{\begin{matrix}{P_{I} = {{- I} + {\frac{\lambda_{1}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}( {B_{2}^{r} - B_{1}^{r}} )} - {\frac{\lambda_{1}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}( {B_{2}^{s} - B_{1}^{s}} )} + ɛ_{I}}} \\{= {{- I} + {\frac{\lambda_{1}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}( {{D\; C\; B^{r}} - {D\; C\; B^{s}}} )} + ɛ_{I}}}\end{matrix}{where}} & (12) \\{N_{I} = {{{- \frac{\lambda_{1}^{3}}{\lambda_{2}^{2} - \lambda_{1}^{2}}}N_{1}} + {\frac{\lambda_{1}^{2}\lambda_{2}}{\lambda_{1} - \lambda_{2}}N_{2}}}} & (13)\end{matrix}$

The relation between L1 ionospheric delay and Total Electron Content is:

$\begin{matrix}{I = {40.3\frac{T\; E\; C}{f_{1}^{2}}}} & (14)\end{matrix}$

where 1 TEC Unit is equivalent to 16.2 cm delay at GPS L1 frequency.

Satellite biases b₁ ^(s) and b₂ ^(s) and receiver biases b₁ ^(r) and b₂^(r) in the ionospheric phase observations cannot be separated from theundifferenced ambiguities, and are thus absorbed in the undifferencedionospheric ambiguities. These biases will be ignored in the followingdescription.

The terms B₂ ^(s)−B₁ ^(s) and B₂ ^(r)−B₁ ^(r) in Eq (6) are theso-called satellite Differential Code Biases and receiver DifferentialCode Biases (DCBs), respectively. The DCB values for GPS and GLONASSsatellites are stable and can be downloaded (e.g., monthly or daily)from the International GNSS Service (IGS) or computed from a globalnetwork of receivers and added to the ionospheric pseudorangeobservation. Further description of the satellite DCBs is given in theBernese Software from the Astronomical Institute, University of Bern,Version 5.0, January 2007, Chapter 13, and in Provisional Applicationfor U.S. Pat. No. 61/337,980 filed 14 Feb. 2010 (TNL A-2633P) andInternational Patent Application PCT/US2011/024733 (TNL A-2633PCT). Thereceiver DCBs need to be estimated in the filter.

United States Patent Application Publication US 2009/0224969 A1 (TNLA-1743US) expresses ionospheric advance across a network area in termsof its Taylor series or other set of orthogonal functions.

Some embodiments of the present invention use a similar model, exceptthat a stochastic ionospheric delay term per satellite per station isadded to account for the inability of describing complex ionosphere witha simple mathematic model:

$\begin{matrix}{{I( {{\Delta\lambda},{\Delta\varphi}} )} = {{{m( {{\Delta\lambda},{\Delta\varphi}} )}( {\sum\limits_{i,{j = 0}}^{\infty}{\alpha_{i,j}{\Delta\lambda}^{i}{\Delta\varphi}^{j}}} )} + \delta}} & (15)\end{matrix}$

where

-   -   Δλ, Δφ are the coordinates of the piercepoints relative to a        reference point,    -   m(Δλ, Δφ) is the so-called mapping function to map ionospheric        slant delay to vertical, and    -   δ is the stochastic delay term which models the small variation        of each satellite-receiver signal path in addition to the Taylor        model.

Since for most purposes the Taylor expansion can be stopped at the firstorder, the terms α_(0,0)=α₀ α_(1,0)=α_(λ) and α_(0,1)=α_(φ) of theTaylor expansion are introduced. The state vector for satellite isub-filter is:

$\begin{matrix}\begin{matrix}{X_{i} = ( {\underset{\underset{{local}\mspace{14mu} {unique}\mspace{14mu} {states}\mspace{14mu} X_{i,u}}{}}{I_{0}^{i},a_{\lambda}^{i},a_{\varphi}^{i},\delta^{1},\ldots \mspace{14mu},\delta^{n},N_{1}^{1},\ldots \mspace{14mu},N_{I}^{n}},\underset{\underset{{Common}\mspace{14mu} {states}\mspace{14mu} X_{l,c}}{}}{B_{I}^{1},{\ldots \mspace{14mu} B_{I}^{n}}}} )^{T}} \\{= ( {X_{i,u},X_{i,c}} )^{T}}\end{matrix} & (16)\end{matrix}$

where states I_(o) ^(i),a_(λ) ^(i),a_(φ) ^(i), . . . , δ^(n),N_(I) ¹, .. . , N_(I) ^(n) are unique in each satellite sub-filter and receiverDCBs B_(I) ¹, . . . B_(I) ^(n) are common for all satellite sub-filter,so that a federated-filter approach can be used.

With the observation update of epoch k in the per satellite sub-filter,the estimated states {circumflex over (X)}_(i) ^(k) and correspondentcovariance matrix Q_(i) ^(k) for satellite i using a standard Kalmanfilter or factorized UD form of Kalman filter are:

$\begin{matrix}{{\hat{X}}_{i}^{k} = \begin{bmatrix}{\hat{X}}_{i,u}^{k} & {\hat{X}}_{i,c}^{k}\end{bmatrix}^{T}} & (17) \\\begin{matrix}{Q_{i}^{k} = \begin{bmatrix}Q_{i,u}^{k} & Q_{i,{uc}}^{k} \\Q_{i,{cu}}^{k} & Q_{i,c}^{k}\end{bmatrix}} \\{= {{\begin{bmatrix}U_{i,u}^{k} & U^{k_{i,{uc}}} \\0 & U_{i,c}^{k}\end{bmatrix}\begin{bmatrix}D_{i,u}^{k} & 0 \\0 & D_{i,c}^{k}\end{bmatrix}}\begin{bmatrix}U_{i,u}^{k} & U_{i,{uc}}^{k} \\0 & U_{i,c}^{k}\end{bmatrix}}^{T}}\end{matrix} & (18)\end{matrix}$

With the federated-filter approach, the receiver DCBs B_(I) ¹, . . .B_(I) ^(n) (common states) and correspondent covariance matrix estimatedin each filter {circumflex over (X)}_(i,c) ^(k), Q_(i,c) ^(k) i=1,2, . .. m are fed into a central fusion master filter.

If a factorized UD form Kalman filter is used as the central fusionmaster filter:

Q_(i,c) ^(k)=U_(i,c) ^(k)D_(i,c) ^(k)U_(i,c) ^(k) ^(T)   (19)

The estimates {circumflex over (X)}_(i,c) ^(k) from single satellite subfilter are first transformed to an uncorrelated vector Z_(i,c) ^(k) withthe inverse of the U matrix:

Z _(i,c) ^(k) =U _(i,c) ^(k) ⁻¹ {circumflex over (X)} _(i,c) ^(k)   (20)

and its corresponding variance is the diagonal matrix D_(i,c) ^(k).

The states in the central fusion master filter (which are the receiverDCBs B_(I) ¹, . . . B_(I) ^(n)) are reset to 0 value and infinitevariance before any observation update at each epoch process; then thede-correlated observations from each satellite sub-filter (eq. 14) areput into a UD filter one by one, from all single-station filters.Storing the covariance matrix in UD form U_(c,k) ^(M), D_(c,k) ^(M)internally, the covariance matrix is computed by:

Q_(c,k) ^(M)=U_(c,k) ^(M)D_(c,k) ^(M)U_(c,k) ^(M) ^(T)   (21)

The central fusion master filter provides global optimized estimates forthe receiver DCBs B_(I) ¹, . . . B_(I) ^(n) which are fed back to eachsatellite's ionospheric sub-filter which then updates the local uniquestates. This assures that the estimates and covariances for the receiverDCBs B_(I) ¹, . . . B_(I) ^(n) in each satellite's ionosphericsub-filter will be the same as in the central fusion master filter.

The local unique states are determined by:

{circumflex over (X)} _(i,u) ^(k) ^(M) ={circumflex over (X)} _(i,u)^(K) −U _(i,uc) ^(k) U _(i,c) ^(k) ⁻¹ ({circumflex over (X)} _(c,k) ^(M)−{circumflex over (X)} _(i,u) ^(K))   (22)

and the covariance matrix for the satellite sub-filter i is:

$\begin{matrix}\begin{matrix}{Q_{i}^{k^{M}} = \begin{bmatrix}Q_{i,u}^{k^{M}} & Q_{i,{uc}}^{k^{M}} \\Q_{i,{cu}}^{k^{M}} & Q_{c,k}^{M}\end{bmatrix}} \\{= {{\begin{bmatrix}U_{i,u}^{k} & U_{i,{uc}}^{k^{M}} \\0 & U_{c,k}^{M}\end{bmatrix}\begin{bmatrix}D_{i,u}^{k} & 0 \\0 & D_{c,k}^{M}\end{bmatrix}}\begin{bmatrix}U_{i,u}^{k} & U_{i,{uc}}^{k^{M}} \\0 & U_{c,k}^{M}\end{bmatrix}}^{T}}\end{matrix} & (23)\end{matrix}$

where

U_(i,uc) ^(k) ^(M) =U_(i,uc) ^(k)U_(i,c) ^(k)U_(c,k) ^(M)   (24)

Eq. (22)-(24) give the final estimates of the local states withcorresponding covariance matrix in UD form. Deriving from eq.(23),

Q_(c,k) ^(M)=U_(c,k) ^(M)D_(i,c) ^(M)U_(c,k) ^(M) ^(T)   (25)

and the covariance for the local unique states is:

$\begin{matrix}{Q_{i,u}^{M} = {{U_{i,u}^{k}D_{i,u}^{k}U_{i,u}^{k}} + {U_{i,{uc}}^{k}D_{c,k}^{M}U_{i,{uc}}^{k^{M^{T}}}}}} & (26)\end{matrix}$

and the covariance between local unique states and common states is:

$\begin{matrix}{Q_{i,{cu}}^{k^{M}} = {U_{c,k}^{M}D_{i,c}^{k^{M}}{U_{i,{uc}}^{k^{M^{T}}}.}}} & (27)\end{matrix}$

PART 3: DERIVING IPWV FROM ZTD

Some embodiments in accordance with the invention use the modifiedHopfield tropospheric model and the Niell (wet and dry) mappingfunction. Other tropospheric models known in the art may be used, suchas the Saastamoinen model, or the Global Pressure and Temperature Modelusing the Global Mapping Function developed by Bohem et al. See, forexample: McCarthy et al., IERS Conventions (2003), IERS Technical Note32, Verlag des Bundesamts für Kartographie and Geodäsie, Frankfurt amMain, 2004, Boehm et al., Troposphere mapping functions for GPS and verylong baseline interferometry from European Centre for Medium-RangeWeather Forecasts operational analysis data, Journal of GeophysicalResearch 111 B02406 DOI:10.1029/2005JB003629 (2006); Boehm et al., TheGlobal Mapping Function (GMF): A new empirical mapping function based ondata from numerical weather model data, Geophysical Research Letters 33L07304 DOI:10.129/2005GL025546 (2006); Boehm et al., Short note: Aglobal model of pressure and temperature for geodetic applications.Journal of Geodesy, Vol. 81, No. 10, pp. 679-683 (2007); and the GlobalPressure and Temperature model published on the IERS ftp siteftp://tai.bipm.org/iers/convupdt/chapter9/gpt.f.

Mapping functions are used to map the zenith tropospheric delay (ZTD) ata receiver to a specific satellite elevation angle. Different mappingfunctions are used for the dry (hydrostatic) component and the wetcomponent:

TD(z)=ƒ_(h)(z)·ZHD+ƒ_(w)(z)·ZWD   (28)

where

TD is the tropospheric Total Delay,

ZHD is the Zenith Hydrostatic Delay,

ZWD is the Zenith Wet Delay,

z is the zenith angle,

ƒ_(h) is the hydrostatic mapping function, and

η_(w) is the wet mapping function.

The Zenith Total Tropospheric delay ZTD represents the total delay,which includes a hydrostatic (often called “dry”) and wet part:

ZTD=ZHD+ZWD   (29)

where ZHD is the zenith hydrostatic delay, and

ZWD is the zenith wet delay.

The Zenith hydrostatic delay ZHD can be accurately computed (tomm-level) if surface pressure is known:

ZHD=0.0022768P/(1−0.00266 cos(2*lat)−0.00028H)   (30)

See Mendes, Modeling the neutral-atmosphere propagation delay inradiometric space techniques. Ph.D. dissertation, Department ofGeomatics Engineering Technical Report No. 199, University of NewBrunswick, Fredericton, New Brunswick, Canada, 353 pp. (1999), availableat http://gge.unb.ca/Pubs/TR199.pdf. Saastamoinen formulas can be usedto evaluate the delay for given pressure. See McCarthy et al., IERSConventions (2003), IERS Technical Note 32, Verlag des Bundesamts fürKartographie and Geodäsie, Frankfurt am Main, 2004, Boehm et al.,Troposphere mapping functions for GPS and very long baselineinterferometry from European Centre for Medium-Range Weather Forecastsoperational analysis data, Journal of Geophysical Research 111 B02406DOI:10.1029/2005JB003629 (2006); Saastamoinen, Atmospheric Correctionfor the Troposphere and Stratosphere in RadioRanging of Satellites,Geophysical Monograph 15, Henriksen (ed), pp. 247-251 (1972); and Daviset al., Geodesy by Radio Interferometry: Effects of AtmosphericModelling Errors on Estimates of Baseline Length, Radio Science, 20, No.6, pp. 1593-1607 (1985).

The additional error in ZHD due to uncertainties in pressure estimationis about 2 mm/mbar for sites not too far from the mean sea level (H=100m and lat=45 deg leads to 0.0023 m/mbar). Pressure values for a givensite can be obtained from various sources, such as local (nearby)measurement, interpolated surface temperature from nearby sites, and/ornumerical weather models.

The Zenith wet delay ZWD can be obtained from GPS-derived ZTD, if ZTD isaccurately available:

ZWD=ZTD−ZHD   (31)

The value of ZWD is highly dependent on the amount of water present inthe atmosphere:

ZWD=10e−6*Rw* (K2+K3/Tm)*IWV   (32)

where IWV is the Integrated Water Vapor—the total mass of water vapor ina 1 m² column extending from the surface to the top of the atmosphere.

Integrated precipitable water vapor (IPWV) is closely related to IWV,where:

IPWV=IWV/ρ_(H) ₂ _(O)   (33)

This means that 1 kg /m² of IWV represents around 1 mm of IPWV. IPWV canbe derived from GNSS observations:

IPWV=(ZTD−ZHD)/[10e−6* Rw*(K2+K3/Tm)*ρ_(H) ₂ _(O)]  (34)

where ZTD is estimated from GNSS measurements,

ZHD is estimated from surface pressure,

Rw, K2 , K3 and ρ_(H) ₂ _(O) are constants,

e is water vapor pressure, and

Tm is mean temperature.

Mean temperature Tm is determined, e.g., from profiles of temperature Tand water vapor e along the atmospheric column, which can come fromradiosonde soundings or numerical weather models.

If temperature T and water vapor e profiles over height z are available,Tm is determined as:

$\begin{matrix}{T_{m} = \frac{\int_{r_{s}}^{r_{a}}{\frac{e}{T}Z_{w}^{- 1}{z}}}{\int_{r_{s}}^{r_{a}}{\frac{e}{T^{2}}Z_{w}^{- 1}{z}}}} & (35)\end{matrix}$

If profiles of temperature T and water vapor pressure e are notavailable , mean temperature Tm can be estimated (less accurately) basedon surface temperatures Ts:

Tm=Ts*[1−bR/(1+1)/g _(m)   (36)

where R is the gas constant for dry air (287.054 J/Kg/K),

g_(m) is the acceleration of gravity computed for given height andlatitude, and

b and l are lapse-rate parameters that can be computed for the locationof interest using a model (e.g., UNB3 m).

Once Tm is known, IPWV can be computed as:

IPWV=(ZTD−ZHD)/[10e−6*Rw*(K2+K3/Tm)*ρ_(H) ₂ _(O)]  (37)

IPWV=ZWD/[10 e−6*Rw*(K2+K3/Tm)*ρ_(H) ₂ _(O)]  (38)

The conversion factor from ZWD to IPWV is strongly dependent on the dayof year and the latitude. Approximate numbers can be estimated in realtime, e.g. every 15 seconds, by a processor using 1 Hz real-time datastreams, with all estimated and measured parameters stored in adatabase. The IPWV information, local surface meteorological data andZWD information can be retrieved from the database as needed. In case ofproblems in real-time data transmission, data can be retrievedautomatically or manually via FTP download from the problematicreference stations, postprocessed, and the results stored in thedatabase. Examples of meteorological data sources include:

Institution Type Format Link NOAA (US) NWM GRIBhttp://nomads.ncdc.noaa.gov/ (Regional, Global) NRCan NWM GRIBhttp://www.weatheroffice.gc.ca/ (CAN) (Regional, grib/ Global) NOAA (US)Radiosonde ASCII http://esrl.noaa.gov/raobs/ (Global) (several)

PART 4: DERIVING TEC/VTEC FROM FIXED AMBIGUITIES

The estimated geometry-free ambiguities from the federated ionosphericfilter and corresponding covariances together with ionospheric-freeambiguities derived from a geometry filter (which may be a federatedgeometry filter as described in United States Patent ApplicationPublication US 2009/0027625 A1 (TNL A-1789US) and widelane ambiguitiesderived from code-carrier filter can be combined to resolve doubledifferenced widelane/narrowlane ambiguities over the network.

With fixed ambiguities, the full accuracy of estimations from carrierphase observations can be reached, i.e. from mm-level to cm-levelaccuracy depending on elevation. Consequently, fixing ambiguities canimprove the estimation of TEC/VTEC and ZTD/Tropospheric slant delay.

Though ambiguities can be expressed in un-differenced,single-differenced, and double differenced form, only thedouble-differenced ambiguities are unique. This means that forun-difference ambiguity, add per satellite bias and per receiver bias,the result double difference ambiguity is unchanged.

The absolute level of ionospheric delay (VTEC) is relatively wellestimated using the federated ionospheric filter approach describedherein. For satellite i and receiver j, for each epoch we can get theestimated float-value ionospheric delay estimates Ĩ_(ν) ^(i,j)=(Ĩ₀^(i),ã_(λ) ^(i),ã_(φ) ^(i),{tilde over (δ)}_(j) ^(i))^(T) from theionospheric filter, and the corresponding covariance matrix Q_(I). Theionospheric slant delay from satellite i to receiver j is then computedas:

$\begin{matrix}\begin{matrix}{{\overset{\sim}{I}}^{i,j} = {\{ {{m( {{\Delta \; \lambda},{\Delta \; \varphi}} )},{{m( {{\Delta \; \lambda},{\Delta \; \varphi}} )} \cdot {\Delta\lambda}},{{m( {{\Delta \; \lambda},{\Delta \; \varphi}} )} \cdot {\Delta\varphi}},1} \} \cdot \begin{Bmatrix}{\overset{\sim}{I}}_{0}^{i} \\{\overset{\sim}{a}}_{\lambda}^{i} \\{\overset{\sim}{a}}_{\varphi}^{i} \\{\overset{\sim}{\delta}}_{j}^{i}\end{Bmatrix}}} \\{= {A \cdot {\overset{\sim}{I}}_{v}^{i,j}}}\end{matrix} & (39)\end{matrix}$

where m(Δλ,Δφ), in(Δλ,Δφ)·Δλ, and m(Δλ,Δφ)·Δφ,1 are the Taylor-expansionmapping functions for the longitudinal and latitudinal gradients. Thecovariance of the ionospheric slant delay is then

σ_(I) ^(i,j) ² =A·Q _(I) ·A ^(T).   (40)

With the (undifferenced) network ambiguities {circumflex over (N)}^(I)^(i,j) fixed for satellite i and receiver j,

{tilde over (L)} _(I) ^(i,j) =L ₁ ^(i,j) −{circumflex over (N)} ₁ ^(i,j)−Ĩ ^(i,j)=Δ^(j) ₁−Δ^(i) ₁+ν₁ ^(i,j),   (41)

where {tilde over (L)}_(I) ^(i,j) is the ambiguity-reduced ionosphericphase observation minus the ionospheric slant delay derived from thefederated ionospheric float filter and serves as input to thefixed-ambiguity ionospheric filter. The terms Δ^(j) _(I) and Δ^(i) _(I)are respectively the receiver biases and satellite biases in theundifferenced ionospheric carrier-phase ambiguities; these cancel if theobservations are double-differenced.

For m satellites and n stations, the number of unknowns is m satellitebiases+n receiver biases. This can be solved either by least squaresadjustment or with a Kalman filter. As the system is rank deficient,either one satellite bias or one receiver bias needs to be constrained;this can be done by introduce a pseudo-observation, for example, byadding an observation Δ_(i) ⁰=0 with a very small variance to define areference satellite bias. After estimating the satellite biases and thestation biases, the improved (corrected) value of the ionospheric slantdelay from the fixed-ambiguity filter for satellite i and receiver j is:

Î ^(i,j) =L ^(i,j) ₁ −{circumflex over (N)} ^(i,j) ₁−{circumflex over(Δ)}_(r) ^(j)+{circumflex over (Δ)}^(i) _(s)   (42)

The values for TEC and VTEC can then be determined in TEC units as

TEC^(i,j) =Î ^(i,j)/0.162   (43)

and

VTEC^(i,j) =m(Δλ,Δφ)·TEC^(i,j).   (44)

PART 5: DERIVING ZTD/TROPOSPHERIC GRADIENT/TROPOSPHERIC SLANT WET DELAYFROM FIXED AMBIGUITIES

The ionospheric-free carrier phase observation and code observation arewritten as:

L _(c) =ρ+c·(t _(r) −t ^(s))+ZTD·ƒ(z)+N _(c)+ν_(c)

P _(c) =ρ+c·(t _(r) −t ^(s))+ZTD·ƒ(z)+ε_(c)   (45)

where t_(r) and t^(s) are receiver clock error and satellite clockerror, respectively,

-   -   ZTD is the Zenith Total Delay,    -   ƒ(z) is the elevation (zenith angle z) dependent mapping        function (in practice, the wet mapping function ƒ_(w) described        in part 3 is taken),    -   N_(c) is the ionospheric free ambiguity, and    -   ν_(c) and ε_(c) are the noise of ionospheric free carrier phase        and code noise.

In case tropospheric gradient needs to be estimated, Eq (45) becomes:

$\begin{matrix}{{{L_{c} = {\rho + {c \cdot ( {t_{r} - t^{s}} )} + {{ZTD} \cdot {f(z)}} + {\frac{\partial f}{\partial z}{{\cos (A)} \cdot x_{N}}} + {\frac{\partial f}{\partial z}{{\sin (A)} \cdot x_{E}}} + N_{c} + v_{c}}}{P_{c} = {\rho + {c \cdot ( {t_{r} - t^{s}} )} + {{ZTD} \cdot {f(z)}} + {\frac{\partial f}{\partial z}{{\cos (A)} \cdot x_{N}}} + {\frac{\partial f}{\partial z}{{\sin (A)} \cdot x_{E}}} + ɛ_{c}}}\mspace{20mu} {{where}\mspace{14mu} \frac{\partial f}{\partial z}}}\;} & (46)\end{matrix}$

is the partial derivative of the mapping function with respect to zenithangle,

A is the azimuth of the direction of station-satellite,

x_(N) and x_(E) are the tropospheric gradients of north and eastdirection, respectively

With m stations and n satellites, the state vector of a geometry filteris

X=[ZTD₁,ZTD₂, . . . ZTD_(m),t_(r1),t_(r2), . . . t_(rm),t^(s1),t^(s2), .. . t^(sn),N_(c1) ¹,N_(c2) ¹, . . . N_(cm) ^(n)]^(T)   (47)

If tropospheric gradient is also estimated, the state vector is

X=[ZTD₁,ZTD₂, . . . ZTD_(m),x_(N1),x_(N2), . . . x_(Nm),x_(E1),x_(E2), .. . x_(Em),t_(r1),t_(r2), . . . t_(rm),t^(s1),t^(s2), . . .t^(sm),N_(c1) ¹,N_(c2) ¹, . . . N_(cm) ^(n)]^(T)   (48)

With an optimal geometry filter or a federated geometry filter, ZTD andTropospheric gradient can be estimated.

With a similar approach as deriving TEC/VTEC from fixed ambiguities,using the network fixed ambiguities, the accuracy of ZTD and slant totaldelay estimation can be improved.

With the network fixed ambiguities {circumflex over (N)}_(c)=N_(c)+Δ^(r)_(c)−Δ^(s) _(c) where Δ^(r) _(c) and Δ^(s) _(c) are respectively thereceiver and satellite dependent biases in the ionosphere-freeundifferenced ambiguities, the ambiguity-reduced, ionosphere-freecarrier phase observation becomes:

$\begin{matrix}\begin{matrix}{{\overset{\Cup}{L}}_{c} = {L_{c} - {\hat{N}}_{c}}} \\{= {\rho + {c \cdot ( {t_{r} - t^{s}} )} + {{ZTD} \cdot {f(z)}} + \Delta_{c}^{s} - \Delta_{c}^{r} + ɛ_{c}}} \\{= {\rho + {c \cdot ( {( {t_{r} - \Delta_{c}^{r}} ) - ( {t^{s} - \Delta_{c}^{s}} )} )} + {{ZTD} \cdot {f(z)}} + ɛ_{c}}} \\{= {\rho + {c \cdot ( ( {{\overset{\Cup}{t}}_{r} - {\overset{\Cup}{t}}^{s}} ) )} + {{ZTD} \cdot {f(z)}} + {ɛ_{c}.}}}\end{matrix} & (49)\end{matrix}$

The terms Δ^(r) _(c) and Δ^(s) _(c) are absorbed by the new satelliteand receiver clock error terms {hacek over (t)}^(s) and {hacek over(t)}_(r). With m satellites and n receivers, the number of unknowns ism+2n , which includes m satellite clock error values plus n receiverclock error values and n ZTD values. Similar to deriving TEC/VTEC fromfixed ambiguities, the system is also rank-deficient and thus either onesatellite clock error or one receiver clock error needs to beconstrained. An improved (corrected) ZTD value for each station can beestimated using least-squares adjustment or a Kalman filter. Optionally,the ZTD values derived from the geometry filter can be used aspseudo-observations.

In case tropospheric gradient needs to be estimated, Eq (49) becomes:

$\begin{matrix}\begin{matrix}{{\overset{\Cup}{L}}_{c} = {L_{c} - {\hat{N}}_{c}}} \\{= {\rho + {c \cdot ( ( {{\overset{\Cup}{t}}_{r} - {\overset{\Cup}{t}}^{s}} ) )} + {{ZTD} \cdot {f(z)}} +}} \\{{{\frac{\partial f}{\partial z}{{\cos (A)} \cdot x_{N}}} + {\frac{\partial f}{\partial z}{{\sin (A)} \cdot x_{E}}} + {v_{c}.}}}\end{matrix} & (50)\end{matrix}$

With m satellites and n receivers, the number of unknowns is m+4n, whichincludes m satellite clock error values plus n receiver clock errorvalues, n ZTD values and 2n tropospheric gradients.

Beside the improved ZTD/tropospheric gradient estimates, slanttropospheric total delay can be derived:

S{circumflex over (T)}D={hacek over (L)}−c·(({hacek over ({circumflexover (t)} ^(r) −{hacek over ({circumflex over (t)} _(s)))   (51)

Where {hacek over ({circumflex over (t)}^(r) and {hacek over({circumflex over (t)}_(s) are the estimates of {hacek over (t)}_(r) and{hacek over (t)}^(s).

With surface meteorogical data, Slant wet delay can be derived:

SŴD=S{circumflex over (T)}D−ƒ_(h)(z)·ZHD−  (52)

Real-time integrated precipitable water vapor (IPWV) can then bedetermined by using estimated ZTD and radiosonde data or surfacemeteorological sensor data as described in part 3.

PART 6: ILLUSTRATIVE EXAMPLES

FIG. 7 is a flow chart of a real-time process 700 for estimatingatmospheric parameter values, using a federated ionospheric filteringmethod in accordance with some embodiments of the invention. Data 705,710, . . . 715 representing measurements of signals from GNSS satellitescollected by receivers of a reference station network distributed over aregion are optionally prepared at 720 (e.g., synchronized by epoch) andmade available as a GNSS data set 725 for processing by epoch. A set ofsatellite DCBs 730, one per satellite, is obtained, e.g., downloadedfrom the International GNSS Service (IGS) or computed from observationsof a global network of receivers such as that shown in FIG. 600 andadded to the ionospheric pseudorange observations.

A federated ionospheric float filter process 735 is applied toionospheric combinations of the GNSS data 725 as adjusted by thesatellite DCBs 730 to estimate a set of receiver DCBs 740, one perstation, and a set of total electron content (TEC) values 745, one perper station per satellite. Process 735 also estimates float values ofionospheric-filter ambiguities 750.

A geometry filter process 755 is applied to the GNSS data to estimatevalues for tropospheric parameters comprising at least one of: (i) a setof values 760 for zenith total delay (ZTD), one per station; (ii) a setof values 760 for zenith total delay (ZTD), one per station, and a setof values 765 for tropospheric gradients per station, one per station ineach of two orthogonal directions (e.g., East and North, or longitudeand latitude), and (iii) a set of values 770 for slant total delay, oneper station per satellite. Process 755 also estimates float values ofgeometry-filter ambiguities 775.

Meteorological data 780 for locations within the region, e.g.,meteorological data from measurements and/or models and optionally fromradiosonde soundings, are obtained from external sources. An atmosphericwatch process 785 determines values over the network region forenvironmental parameters 790. At least one of (1) integratedprecipitable water vapor (IPWV) per station 792 and (2) troposphericslant wet delay 794 is determined from the estimated values fortropospheric parameters and the meteorological data 780. The atmosphericwatch process or another process optionally determines a set of values796 for vertical total electron content (VTEC), one per station, fromthe values of total electron content (TEC) per station per satellite745.

FIG. 8 schematically illustrates an architecture using a federatedionospheric filter in accordance with some embodiments of the invention.A federated ionospheric float filter 835 is applied to ionosphericcombinations of the GNSS data 725 as adjusted by the satellite DCBs 730to estimate a set of receiver DCBs 740, one per station, and a set oftotal electron content (TEC) values 745, one per station per satellite.Federated filter 835 also estimates float values of ionospheric-filterambiguities 750.

A geometry filter 855 is applied to the GNSS data to estimate values fortropospheric parameters comprising at least one of: (i) a set of values760 for zenith total delay (ZTD), one per station; (ii) a set of values760 for zenith total delay (ZTD), one per station, and a set of values765 for tropospheric gradients per station, one per station in each oftwo orthogonal directions (e.g., East and North, or longitude andlatitude), and (iii) a set of values 770 for slant total delay, one perstation per satellite. Filter 855 also estimates float values ofgeometry-filter ambiguities 775.

An atmospheric watch module 885 determines values over the networkregion for environmental parameters 790. Module 885 determines values ofat least one of (1) integrated precipitable water vapor (IPWV) 792 and(2) tropospheric slant wet delay 794, from the estimated values fortropospheric parameters and from the meteorological data 780. Optionallythe atmospheric watch module or another module determines a set ofvalues 796 for vertical total electron content (VTEC), one per station,from the values of total electron content (TEC) per station persatellite 745.

FIG. 9 is a flow chart of federated ionospheric filter process 735 inaccordance with some embodiments of the invention. The satellite DCBs730 and the GNSS data set 725 are used in a process 905 which estimatesvalues 910 for parameters of an ionospheric model and values 915 of astochastic ionospheric delay term per satellite per station. Theestimated values 910 and 915 are used in a process 920 to determinevalues 745 of total electron content (TEC) per station per satellite.

FIG. 10 is a schematic diagram of an embodiment of federated ionosphericfilter 835 in accordance with some embodiments of the invention. Thesatellite DCBs 730 and the GNSS data set 725 are used in an estimator1005 determines values 910 for parameters of an ionospheric model andvalues 915 of a stochastic ionospheric delay term per satellite perstation. The estimated values 910 and 915 are passed to a module 1020which determines values 745 of total electron content (TEC) per stationper satellite.

FIG. 11 is a flow chart of federated ionospheric filter process 735 inaccordance with some embodiments of the invention. The process begins anepoch of processing at 1005. The GNSS data set 725 and satellite DCBs730 are applied at 1110 to an ionospheric sub-filter process persatellite 1110. Each sub-filter process 1110 estimates values 1115 for“local” states unique to its assigned satellite and values 1120 forstates common to all ionospheric sub-filters. As shown in Eq. (16), the“local” state values 1115 comprise parameters (I₀ ^(i),a_(λ) ^(i),a_(φ)^(i)) characterizing the ionosphere across the region, a stochasticionospheric delay term per satellite per station (δ¹, . . . , δ^(n)),and float-value ambiguities per satellite per station (N_(I) ¹, . . .N_(I) ^(n)). The common state values 1120 comprising a receiver DCB perstation (B_(I) ¹, . . . B_(I) ^(n)), are fed one by one to a centralmaster fusion filter process 1135 which estimates updated values 1145for the common states. An ionospheric sub-filter process 1150 (e.g.,using the same sub-filters a process 1110) uses the updated common statevalues 1145 to estimate updated values 1155 for the “local” states.Feeding back the updated common state values 1145 to the sub-filtersdeals with cross-correlation between the local states and the commonstates.

FIG. 12 is a schematic diagram chart of federated ionospheric filter 835in accordance with some embodiments of the invention. Ionosphericsub-filters 1210, one per satellite, estimate values 1115 for “local”states unique to the each sub-filter's assigned satellite and values1120 for states common to all ionospheric sub-filters. As shown in Eq.(16), the “local” state values 1115 comprise parameters (I_(o)^(i),a_(λ) ^(i),a_(φ) ^(i)) characterizing the ionosphere across theregion, a stochastic ionospheric delay term per satellite per station(δ¹, . . . , δ^(n)), and float-value ambiguities per satellite perstation (N_(I) ¹, . . . , N_(I) ^(n)). The common state values 1120comprising a receiver DCB per station (B_(I) ¹, . . . B_(I) ^(n)), arefed one by one to a central master fusion filter 1235 which estimatesupdated values 1145 for the common states. The ionospheric sub-filters1210 use the updated common state values 1145 to estimate updated values1155 for the “local” states. Feeding back the updated common statevalues 1145 to the sub-filters deals with cross-correlation between thelocal states and the common states.

FIG. 13 is a flow chart of atmospheric watch process 785 in accordancewith some embodiments of the invention. Zenith dry delay (ZHD) perstation is computed at 1305 from meteorological data 1310, e.g., fromdata representing surface meteorological conditions near the station, asin Eq. (30). At 1315, Zenith wet delay (ZWD) per station 1330 andtropospheric slant wet delay per station per satellite 794 are computedfrom zenith dry delay (ZHD) 1315, zenith total delay (ZTD) 760, slanttotal delay 770 and the hydrostatic (dry) mapping function 1325, as inEq. (28). The hydrostatic mapping function is a dimensionless factorwhich describes the elevation angle dependence of the hydrostatic pathdelay and relates the line of sight delay to the zenith delay. At 1335,mean temperature 1345 is computed from radiosonde data 1340 and surfacemeteorological data 1310 (and/or data from numerical weather models), asin Eq. (35) or Eq. (36). Optionally, at 1350, values 792 for integratedprecipitable water vapor (IPWV) are computed from the mean temperature1345 and zenith wet delay (ZWD) 1330 as in Eq. (34). Optionally, values796 of vertical total electron content per station (VTEC) are computedat 1355 from the values 745 of total electron content (TEC) per stationper satellite. As described above, the relationship betweengeometry-dependent TEC and geometry-independent VTEC is given by

VTEC=TEC/m _(mapping)(ζ)=TEC·cos(ζ)   (53)

where m_(mapping)(ζ)=1/cos(ζ) is the slant path mapping function, and

-   -   ζ is the angle between the signal ray path from satellite to        receiver and the line perpendicular to the ionospheric sphere at        the piercepoint of the signal ray path through the ionospheric        sphere.        The environmental parameter estimates 790 (at least one of VTEC        per station 796, IPWV 792, and slant wet delay per station per        satellite 794) are stored in memory 1365, from which they can be        communicated to one or more computer systems, e.g, computer        system 1370, for display and/or further processing such as        weather forecasting.

FIG. 14 is a schematic diagram of atmospheric watch module 885 inaccordance with some embodiments of the invention. A zenith dry delaymodule 1405 computes zenith dry delay (ZHD) from meteorological data1310. A zenith wet delay & slant wet delay module 1420 computes zenithwet delay (ZWD) per station 1330 and tropospheric slant wet delay perstation per satellite 794. A mean temperature module 1435 computes meantemperature 1345. Optionally, an IPWV module 1450 computes values 792for integrated precipitable water vapor (IPWV). Optionally, a VTECmodule 1455 computes values 796 of vertical total electron content perstation (VTEC). Data processor & delay hardware includes memory 1365 tostore the environmental parameter estimates 790, from which they can becommunicated to one or more computer systems such as computer system1370 for display and/or further processing such as weather forecasting.

FIG. 15 is a flow chart of a process 1500 for estimating environmentalparameters using fixed ambiguities in accordance with some embodimentsof the invention. Data 1505, 1510, . . . 1515 representing measurementsof signals from GNSS satellites collected by receivers of a referencestation network distributed over a region are optionally prepared at1520 (e.g., synchronized by epoch) and made available as a GNSS data set1525 for processing by epoch. A set of satellite DCBs 1530, one persatellite, is obtained, e.g., downloaded from the International GNSSService (IGS) or computed from observations of a global network ofreceivers such as that shown in FIG. 600 and added to the pseudorangeobservations.

At 1535 a recursive filtering process (using e.g., a Least Squaresadaptive filter or Kalman filter) is applied to the GNSS data 1525 asadjusted by the satellite DCBs 1530 to estimate values at each epoch fora set of state variables 1540. State variables 1540 include at least oneof: (i) total electron count (TEC) per station per satellite, (ii) azenith total delay (ZTD) per station, (iii) a zenith total delay (ZTD)per station and a set of tropospheric gradients per station, and (iv) aslant total delay per station per satellite. The process also estimatesfloat values for ambiguities and associated statistical information(variance/covariance values).

At 1545, the float-value ambiguity estimates are “fixed” to form a setof “fixed” network ambiguities 1550. At 1555, the “fixed” networkambiguities 1550 are used to estimate corrected values 1560 for thestate variables. Optionally, “fixed” network ambiguities 1550 can be fedback to Recursive Filter 1635. The resulting values for atmosphericparameters (at least one of (i) corrected total electron count (TEC) perstation per satellite 1565, (ii) corrected zenith total delay (ZTD) perstation 1570, (iii) corrected zenith total delay (ZTD) per station 1570and corrected tropospheric gradients per station 1575, and (iv)corrected slant total delay per station per satellite 1580) andmeteorological data 1585 from external sources, are supplied to anatmospheric watch process, e.g., as at 785 in FIG. 7. The atmosphericwatch process determines environmental parameter values 1590 over theregion for at least one of (1) integrated precipitable water vapor(IPWV) 1592 and (2) tropospheric slant wet delay 1594. Themeteorological data 1585 for locations within the region comprise, e.g.,meteorological data from measurements and/or models and optionally fromradiosonde soundings. The atmospheric watch process or another processdetermines values 1596 of vertical total electron content (VTEC) perstation from the values of total electron content per station persatellite (TEC) 1565.

FIG. 16 is a schematic diagram of a system for estimating environmentalparameter values using fixed ambiguities in accordance with someembodiments of the invention. An optional data preparation module 1620prepares (e.g., synchronizes) data 1505, 1510, . . . 1515 representingmeasurements of signals from GNSS satellites collected by receivers of areference station network distributed over a region and makes theprepared data available as a GNSS data set 1525 for processing by epoch.A recursive filter 1635 (e.g., one or more Least Squares adaptivefilters or Kalman filters) is applied to the GNSS data 1525 as adjustedby satellite DCBs 1530 to estimate values at each epoch for a set ofstate variables 1540. State variables 1540 include at least one of (i) atotal electron count (TEC) per station per satellite, (ii) a zenithtotal delay (ZTD) per station, (iii) a zenith total delay (ZTD) perstation and a set of tropospheric gradients per station, and (iv) aslant total delay per station per satellite. Filter 1635 also estimatesfloat values for ambiguities and associated statistical information(variance/covariance values).

An ambiguities “fixing” module 1645 “fixes” the float-value ambiguityestimates to form a set of “fixed” network ambiguities 1550. Astate-variable estimates corrector 1655 uses the “fixed” networkambiguities 1550 to estimate corrected values 1560 for the statevariables. Optionally, atmospheric watch module 885 uses the resultingcorrected values 1560 for the state variables (at least one of: (i)corrected total electron count (TEC) per station per satellite 1565,(ii) corrected zenith total delay (ZTD) per station 1570, (iii)corrected zenith total delay (ZTD) per station 1570 and correctedtropospheric gradients per station 1575, and (iv) corrected slant totaldelay per station per satellite 1580), and meteorological data 1585 fromexternal sources, to determine environmental parameter values 1590 overthe region for at least one of (1) integrated precipitable water vapor(IPWV) 1592 and (2) tropospheric slant wet delay 1594. Themeteorological data 1585 for locations within the region comprise, e.g.,meteorological data from measurements and/or models and optionally fromradiosonde soundings. Optionally, the atmospheric watch process oranother process determines values 1596 of vertical total electroncontent (VTEC) per station from the values of total electron content perstation per satellite (TEC) 1565.

FIG. 17 is a flow chart of a process 1700 for estimating environmentalparameters using fixed ambiguities in accordance with some embodimentsof the invention. Data 1705, 1710, . . . 1715 representing measurementsof signals from GNSS satellites collected by receivers of a referencestation network distributed over a region are optionally prepared at1720 (e.g., synchronized by epoch) and made available as a GNSS data set1725 for processing by epoch. A set of satellite DCBs 1730, one persatellite, is obtained, e.g., downloaded from the International GNSSService (IGS) or computed from observations of a global network ofreceivers such as that shown in FIG. 600 and added to the pseudorangeobservations.

At 1735, recursive filtering is performed using factorized filterprocesses applied to the GNSS data 1725 to estimate values for a set ofstate variables at each epoch. An ionospheric process 1740 uses anionospheric (geometry-free) observation combination, adjusted by thesatellite DCBs 1730, to estimate values 1742 for a receiver differentialcode bias (DCB) per station, optionally values 1744 for total electroncount per station per satellite, and float values 1746 for theionospheric-combination ambiguities, along with associated statisticalinformation. A geometry-filter process 1750 uses an ionospheric-free(geometric) observation combination to estimate values for atmosphericparameters (comprising at least one of zenith total delay (ZTD) perstation 1752, zenith total delay (ZTD) per station and values 1754 fortropospheric gradient per station, values 1756 for slant total delay perstation), and float values 1750 for the geometric-combinationambiguities, along with associated statistical information. Acode/carrier filter process 1760 uses a code/carrier observationcombination to estimate values 1762 for the code/carrier parameters andfloat values 1764 for the code/carrier-combination ambiguities, alongwith associated statistical information. At 1765 the estimates and theassociated statistical information are combined to form a combined set1770 of estimates with associated statistical information (e.g., avariance/covariance matrix).

Factorized filtering can be performed as described in U.S. Pat. No.7,432,853 (TNL A-1403US), except that here the GNSS data applied to theionospheric filter process is adjusted by the satellite DCBs. In someembodiments, the ionospheric filtering of U.S. Pat. No. 7,432,853 isreplaced with the federated ionospheric filtering process describedabove. In some embodiments, the geometry filtering process of U.S. Pat.No. 7,432,853 is replaced with the federated geometry filter process ofUnited States Patent Application Publication US 2009/0027625 A1 (TNLA-1789US).

The combined array 1770 of ambiguity estimates for all carrier phaseobservations and associated statistical information are used at 1775to“fix” the ambiguity values. The resulting “fixed” network ambiguities1776 are used at 1780 to estimate corrected values for the environmentalstate variables. The resulting corrected values for the environmentalstate variables (at least one of (i) corrected total electron count(TEC) per station per satellite 1782, (ii) corrected values 1784 ofzenith total delay (ZTD) per station, (iii) corrected values 1784 ofzenith total delay (ZTD) per station and corrected values 1786 of a setof tropospheric gradients per station, and (iv) corrected values 1788 ofslant total delay per station per satellite), together withmeteorological data 1785 from external sources, are optionally suppliedto an atmospheric watch process, e.g., process 785 of FIG. 7. Optionalprocess 785 uses the corrected values to determine values 1792 for atleast one of (1) integrated precipitable water vapor (IPWV), (2) values1794 for tropospheric slant wet delay, and (3) values 1796 for verticaltotal electron count (VTEC). The meteorological data 1785 comprise,e.g., meteorological data from measurements at locations within theregion and/or models and optionally from radiosonde soundings.

FIG. 18 is a schematic diagram of a system for estimating environmentalparameters using fixed ambiguities in accordance with some embodimentsof the invention. A data preparation module 1820 prepares (e.g.,synchronizes) data 1705, 1710, . . . 1715 representing measurements ofsignals from GNSS satellites collected by receivers of a referencestation network distributed over a region and makes these available as aGNSS data set 1725 for processing by epoch. A factorized filter 1835 isuses the GNSS data 1725 to estimate values for a set of state variablesat each epoch. Factorized filter 1835 includes one or more ionosphericfilters 1840 which use an ionospheric (geometry-free) observationcombination, adjusted by the satellite DCBs 1730, to estimate values1742 for a receiver differential code bias (DCB) per station, optionallyvalues 1744 for total electron count per station per satellite, andfloat values 1746 for the ionospheric-combination ambiguities, alongwith associated statistical information. Factorized filter 1835 includesa geometry-filter 1850 which uses an ionospheric-free (geometric)observation combination to estimate values for atmospheric parameters(comprising at least one of (i) zenith total delay (ZTD) per station1752, (ii) zenith total delay (ZTD) per station 1752 and troposphericgradient per station 1754, and (iii) slant total delay per station1754), and float values 1750 for the geometric-combination ambiguities,along with associated statistical information. Factorized filter 1835includes code/carrier filters 1860 which use a code/carrier observationcombination to estimate values 1762 for the code/carrier parameters andfloat values 1764 for the code/carrier-combination ambiguities, alongwith associated statistical information. A combiner 1865 combines theestimates and the associated statistical information into a combined set1770 of estimates with associated statistical information (e.g., avariance/covariance matrix).

Factorized filter 1835 can be implemented as described in U.S. Pat. No.7,432,853 (TNL A-1403US), except that here the GNSS data (theionospheric combination) applied to the ionospheric filter(s) 1840 isadjusted by the satellite DCBs. In some embodiments, the ionosphericfilters 1840 of U.S. Pat. No. 7,432,853 are replaced with the federatedionospheric filter described above. In some embodiments, the geometryfilter of U.S. Pat. No. 7,432,853 is replaced with the federatedgeometry filter of United States Patent Application Publication US2009/0027625 A1 (TNL A-1789US).

An ambiguities “fixing” module 1875 “fixes” the ambiguities of thecombined array 1770. A corrector 1880 uses the “fixed” networkambiguities 1776 to estimate corrected values for the environmentalstate variables. Optionally, atmospheric watch module 885 uses thecorrected environmental state values (at least one of (i) values 1782 oftotal electron count (TEC) per station per satellite, (ii) correctedvalues 1784 of zenith total delay (ZTD) per station, (iii) correctedvalues 1784 of zenith total delay (ZTD) per station and corrected values1786 of tropospheric gradients per station, and (iv) corrected values1788 of slant total delay per station per satellite), together withmeteorological data 1785 from external sources, to determineenvironmental parameter estimates 1790. Optional estimates 1790 includevalues for at least one of (1) integrated precipitable water vapor(IPWV) 1792, (2) tropospheric slant wet delay 1794, and (iv) verticaltotal electron count (VTEC) 1796. The meteorological data 1785 comprise,e.g., meteorological data from measurements at locations within theregion and/or models and optionally from radiosonde soundings.

The terms “fix,” “fixed” and “fixing” as used herein include fixing ofambiguities to integer values using techniques such as rounding,bootstrapping and Lambda search, and also include forming a weightedaverage of integer candidates while not necessarily fixing them tointeger values. The weighted average approach is described in detail inInternational Patent Publications WO 2010/021656 A2, WO 2010/021658 A2,WO 2010/021659 A2, WO 2010/021657 A2, and WO 2010/021660 A2,incorporated herein by reference. Some embodiments “fix” ambiguitiesusing any suitable techniques known in the art, such as Simple rounding,Bootstrapping , Integer Least Squares based on the Lambda Method, orBest Integer Equivariant. See Teunissen, P. J. G, S. Verhagen (2009);GNSS Carrier Phase Ambiguity Resolution: Challenges and Open Problems,in M. G. Sideris (ed.); Observing our Changing Earth, InternationalAssociation of Geodesy Symposia 133, Springer Verlag Berlin-Heidelberg2009; and Verhagen, Sandra (2005): The GNSS integer ambiguities:estimation and validation, Publications on Geodesy 58, Delft, 2005, 194pages, ISBN-13: 978 90 6132 290 0. ISBN-10: 90 6132 290 1.

-   -   N₁ and N₂ are L₁ and L₂ integer ambiguities,    -   ν₁ and ν₂ are L₁ and L₂ phase noise plus multipath, and    -   ε_(l) and ε₂ are L₁ and L₂ code noise plus multipath.

The carrier-phase observations of Eqs. (7) and (8) show the L1 and L2ambiguities as N₁ and N₂, respectively. Eq. (7) shows the mappingbetween N₁ and N₂. For a network of n stations observing m satellites,the m x n ambiguities on each of the carrier frequencies, e.g., on L1and L2, can be expressed in various ways. Linear combinations ofobservations result in corresponding combinations of the ambiguitiesfrom which the uncombined ambiguities can be derived. For example, asingle large recursive filter as in the example of FIG. E can beconfigured to directly estimate the uncombined ambiguities. In theexample of FIG. F using factorized filters, each of the factorizedfilters estimates a different linear combination of the ambiguities,e.g., widelane, narrowlane, and code/carrier combinations. Alternately,the ambiguities can be expressed as baseline-to-baseline combinations,where each baseline corresponds to a respective pair of the networkreceivers. The present invention is intended to include estimation ofany combinations of ambiguities which are equivalent to fixinguncombined ambiguities across the network. Discussion of some of thelinear combinations is found, for example, in Chapter 2 of BerneseSoftware from the Astronomical Institute, University of Bern, Version5.0, January 2007.

PART 7: RECEIVER AND PROCESSING APPARATUS

FIG. 19 is a block diagram of a typical integrated GNSS receiver system1900 with GNSS antenna 1905 and communications antenna 1910. The TrimbleR8 GNSS System is an example of such a system. Receiver system 1900 canserve as a network reference station. Receiver system 1900 includes aGNSS receiver 1915, a computer system 1920 and one or morecommunications links 1925. Computer system 1920 includes one or moreprocessors 1930, one or more data storage elements 1935, program code1940 with instructions for controlling the processor(s) 1930, and userinput/output devices 1945 which may include one or more output devices1950 such as a display or speaker or printer and one or more devices1955 for receiving user input such as a keyboard or touch pad or mouseor microphone.

FIG. 20 is a schematic diagram of a computer system 2020 in accordancewith some embodiments of the invention. Computer system 2020 includesone or more processors 2030, one or more data storage elements 2035,program code 2040 with instructions for controlling the processor(s)2030, and user input/output devices 2045 which may include one or moreoutput devices 2050 such as a display or speaker or printer and one ormore devices 2055 for receiving user input such as a keyboard or touchpad or mouse or microphone. Computer system 2020 may also have one ormore communications links for exchanging data with other systems, e.g.,via internet and/or other channels.

PART 8: GENERAL REMARKS

The inventive concepts can be employed in a wide variety of processesand equipment for real-time operation (no more than 5 seconds). Usingthe federated filter approach and factorized UD filter greatly reducesthe computation load for large GNSS network and is suitable to run withmodern multi-core computer with parallel computation technique(i.e. byusing OpenMP library). Shorter time delay is crutial to life criticalapplication like flood and weather predition near an airport. Someexemplary embodiments will now be described. It will be understood thatthese are intended to illustrate rather than to limit the scope of theinvention.

Those of ordinary skill in the art will realize that the detaileddescription of embodiments of the present invention is illustrative onlyand is not intended to be in any way limiting. Other embodiments of thepresent invention will readily suggest themselves to such skilledpersons having the benefit of this disclosure. For example, while aiono-free phase combination is employed in some examples, those of skillin the art will recognized that many phase combinations are possible andthat a combination other than a iono-free phase combination can produceacceptable if less than optimum results; thus the claims are notintended to be limited to iono-free phase combinations other than whereexpressly called for. Reference is made in detail to implementations ofthe present invention as illustrated in the accompanying drawings. Thesame reference indicators are used throughout the drawings and thefollowing detailed description to refer to the same or like parts.

In the interest of clarity, not all of the routine features of theimplementations described herein are shown and described. It will beappreciated that in the development of any such actual implementation,numerous implementation-specific decisions must be made to achieve thedeveloper's specific goals, such as compliance with application- andbusiness-related constraints, and that these specific goals will varyfrom one implementation to another and from one developer to another.Moreover, it will be appreciated that such a development effort might becomplex and time-consuming, but would nevertheless be a routineundertaking of engineering for those of ordinary skill in the art havingthe benefit of this disclosure.

Global Navigation Satellite Systems (GNSS) include the GlobalPositioning System (GPS), the Glonass system, the proposed Galileosystem, the proposed Compass system, and others. Each GPS satellitetransmits continuously using two radio frequencies in the L-band,referred to as L1 and L2, at respective frequencies of 1575.41 MHz and1227.60 MHz. Some GPS satellites also transmit on a third radiofrequency in the L-band, referred to as L5 at 1176.45 MHz. Each GNSSlikewise has satellites which transmit multiple signals on multiplecarrier frequencies. Embodiments of the present invention are notlimited to any specific GNSS, or to L1 and L2 frequencies. For example,the invention may be used with other combinations of frequencies, suchas GPS L1 and L5, or combinations of Glonass frequencies, combinationsof Gallileo frequencies, or combinations of frequencies of any otherGNSS.

In accordance with embodiments of the present invention, the components,process steps and/or data structures may be implemented using varioustypes of operating systems (OS), computer platforms, firmware, computerprograms, computer languages and/or general-purpose machines. Themethods can be run as a programmed process running on processingcircuitry. The processing circuitry can take the form of numerouscombinations of processors and operating systems, or a stand-alonedevice. The processes can be implemented as instructions executed bysuch hardware, by hardware alone, or by any combination thereof. Thesoftware may be stored on a program storage device readable by amachine. Computational elements, such as filters and banks of filters,can be readily implemented using an object-oriented programming languagesuch that each required filter is instantiated as needed.

Those of skill in the art will recognize that devices of a lessgeneral-purpose nature, such as hardwired devices, field programmablelogic devices (FPLDs), including field programmable gate arrays (FPGAs)and complex programmable logic devices (CPLDs), application specificintegrated circuits (ASICs), or the like, may also be used withoutdeparting from the scope and spirit of the inventive concepts disclosedherein.

In accordance with an embodiment of the present invention, the methodsmay be implemented on a data processing computer such as a personalcomputer, workstation computer, mainframe computer, or high-performanceserver running an OS such as Microsoft® Windows® XP and Windows® 2000,available from Microsoft Corporation of Redmond, Wash., or Solaris®available from Sun Microsystems, Inc. of Santa Clara, Calif., or variousversons of the Unix operating system such as Linux available from anumber of vendors. The methods may also be implemented on amultiple-processor system, or in a computing environment includingvarious peripherals such as input devices, output devices, displays,pointing devices, memories, storage devices, media interfaces fortransferring data to and from the processor(s), and the like. Such acomputer system or computing environment may be networked locally, orover the Internet.

Any of the above-described methods and their embodiments may beimplemented by means of a computer program. The computer program may beloaded on an apparatus, a rover, a reference receiver or a networkstation as described above. Therefore, the invention also relates to acomputer program, which, when carried out on an apparatus, a rover, areference receiver or a network station as described above, carries outany one of the above above-described methods and their embodiments.

The invention also relates to a computer-readable medium or acomputer-program product including the above-mentioned computer program.The computer-readable medium or computer-program product may forinstance be a magnetic tape, an optical memory disk, a magnetic disk, amagneto-optical disk, a CD ROM, a DVD, a CD, a flash memory unit or thelike, wherein the computer program is permanently or temporarily stored.The invention also relates to a computer-readable medium (or to acomputer-program product) having computer-executable instructions forcarrying out any one of the methods of the invention.

The invention also relates to a firmware update adapted to be installedon receivers already in the field, i.e. a computer program which isdelivered to the field as a computer program product. This applies toeach of the above-described methods and apparatuses.

GNSS receivers may include an antenna, configured to received thesignals at the frequencies broadcasted by the satellites, processorunits, one or more accurate clocks (such as crystal oscillators), one ormore computer processing units (CPU), one or more memory units (RAM,ROM, flash memory, or the like), and a display for displaying positioninformation to a user.

Where the terms “receiver”, “filter” and “processing element” are usedherein as units of an apparatus, no restriction is made regarding howdistributed the constituent parts of a unit may be. That is, theconstituent parts of a unit may be distributed in different software orhardware components or devices for bringing about the intended function.Furthermore, the units may be gathered together for performing theirfunctions by means of a combined, single unit. For instance, thereceiver, the filter and the processing element may be combined to forma single unit, to perform the combined functionalities of the units.

The above-mentioned units may be implemented using hardware, software, acombination of hardware and software, pre-programmed ASICs(application-specific integrated circuit), etc. A unit may include acomputer processing unit (CPU), a storage unit, input/output (I/O )units, network connection units, etc.

Although the present invention has been described on the basis ofdetailed examples, the detailed examples only serve to provide theskilled person with a better understanding, and are not intended tolimit the scope of the invention. The scope of the invention is muchrather defined by the appended claims.

PART 9: SUMMARY OF INVENTIVE CONCEPTS

Following is a summary of some of the inventive concepts describedherein

Part 9A: [GNSS Atmospheric Estimation with Federated Ionospheric Filter]

-   -   1. A method of processing GNSS signal data to estimate        environmental parameter values, comprising        -   obtaining GNSS data collected at a plurality of stations            distributed over a region from signals received from GNSS            satellites over multiple epochs,        -   obtaining a satellite differential code bias (DCB) per            satellite, and        -   estimating a receiver differential code bias per station and            a value of total electron content (TEC) per station per            satellite from the GNSS data and the satellite differential            code biases using a federated ionospheric filter.    -   2. The method of 1, further comprising        -   applying a geometric filter to estimate from the GNSS data a            set of values for atmospheric parameters comprising at least            one of: (i) a zenith total delay per station, (ii) a zenith            total delay per station and a set of tropospheric gradients            per station, and (iii) a slant total delay per station per            satellite,        -   obtaining meteorological data for locations within the            region, and        -   determining values over the region for at least one of: (1)            integrated precipitable water vapor (IPWV) and (2)            tropospheric slant wet delay, from the set of values for            atmospheric parameters and the meteorological data.    -   3. The method of 1 or 2, wherein obtaining a satellite        differential code bias per satellite comprises one of (i)        retrieving a satellite differential code bias per satellite from        an external source, and (ii) computing a satellite differential        code bias per satellite from GNSS data collected at a network of        reference stations.    -   4. The method of one of 1-3, wherein estimating a value of total        electron content per station per satellite comprises        -   estimating values for ionospheric model parameters and for a            stochastic’ ionospheric delay term per satellite per station            in the federated ionospheric filter, and        -   determining the value of total electron content per station            per satellite from the estimated values for the ionospheric            model parameters and the stochastic ionospheric delay terms.    -   5. The method of one of 1-4, wherein applying a federated        ionospheric filter comprises:        for each satellite,    -   applying to the GNSS signal data an ionospheric subfilter to        estimate values for local states representing parameters unique        to that satellite and values for common states representing        parameters common to all receivers,        -   -   providing values for the common states and associated                statistical information to a master filter; and            -   preparing updated estimates for the local states when                updated values for the common states are provided by the                master filter; and

        -   applying to the values for the common states and the            associated statistical information a master filter to            estimate updated values for the common states, and to            provide the updated values to the ionospheric subfilters.    -   6. The method of 5, wherein the states unique to the satellite        comprise: a set of parameters characterizing ionosphere across        the region, a stochastic ionospheric delay term per station per        satellite, and an integer ambiguity per station per satellite.    -   7. The method of one of 5-6, wherein the states common to all        satellites comprise a receiver differential code bias per        station.    -   8. The method of 5, further comprising: resetting states in the        master filter to zero with infinite variance before any        observation update at each epoch, and then applying the        decorrelated observations from each subfilter to the master        filter    -   9. The method of one of 1-8, wherein the value of total electron        content per station per satellite comprises a value mapped to        vertical.    -   10. The method of one of 2-9, wherein obtaining meteorological        data comprises obtaining surface meteorological data for        locations within the region.    -   11. The method of one of 2-10, wherein obtaining meteorological        data comprises obtaining radiosonde temperature data for        locations within the region.    -   12. The method of one of 1-11, wherein elapsed time between        obtaining the GNSS data and determining values over the region        for total electron content is no more than about five seconds.    -   13. The method of one of 2-11, wherein elapsed time between        obtaining the GNSS data and determining values over the region        for total electron content and at least one of integrated        precipitable water vapor and tropospheric slant delay is no more        than about five seconds.    -   14. Apparatus for processing GNSS signal data to estimate        environmental parameter values, comprising a federated        ionospheric filter adapted to estimate a receiver differential        code bias per station and a value of total electron content        (TEC) per station per satellite from satellite differential code        biases and the GNSS signal data.    -   15. The apparatus of 14, further comprising:        -   a geometric filter adapted to estimate from the GNSS data a            set of values for atmospheric parameters comprising at least            one of: (i)a zenith total delay per station, (ii) a zenith            total delay per station and a set of tropospheric gradients            per station, and (iii)a slant total delay per station per            satellite, and        -   an atmospheric watch module to determine values over the            region for integrated precipitable water vapor and            tropospheric slant wet delay, from the estimated values for            atmospheric parameters and meteorological data for the            region.    -   16. The apparatus of one of 14-15, further comprising an element        to obtain a code bias per satellite by one of (i) retrieving a        code bias per satellite from an external source, and (ii)        computing a code bias per satellite from GNSS data collected at        a network of GNSS reference stations.    -   17. The apparatus of one of 14-16, comprising:        -   a module to estimate values for ionospheric model parameters            and for a stochastic ionospheric delay term per satellite            per station, and        -   a module to determine the value of total electron content            per station per satellite from the estimated values for the            ionospheric model parameters and the stochastic ionospheric            delay terms.    -   18. The apparatus of one of 14-17, wherein the federated        ionospheric filter comprises:        -   an ionospheric subfilter for each satellite adapted to            -   estimate values for local states' representing                parameters unique to that satellite and values for                common states representing parameters common to all                receivers, from the GNSS signal data and a satellite                differential code bias,            -   provide values for the common states and associated                statistical information to master filter; and            -   prepare updated estimates for the local states when                updated values for the common states are provided by the                master filter; and        -   a master filter adapted to estimate updated values for the            common states from the values for the common states and the            associated statistical information, and to provide the            updated values to the ionospheric subfilters.    -   19. The apparatus of 18, wherein the states unique to the        satellite comprise: a set of parameters characterizing        ionosphere across the region, a stochastic ionospheric delay        term per station per satellite, and an integer ambiguity per        station per satellite.    -   20. The apparatus of one of 18-19, wherein the states common to        all satellites comprise a receiver differential code bias per        station.    -   21. The apparatus of 18, wherein the federated ionospheric        filter is operative to reset states in the master filter to zero        with infinite variance before observation update at each epoch,        and then to apply the decorrelated observations from each        subfilter to the master filter.    -   22. The apparatus of one of 14-21, wherein the value of total        electron content per station per satellite comprises a value        mapped to vertical.    -   23. The apparatus of one of 14-22, wherein the meteorological        watch module is adapted to use surface meteorological data for        locations within the region as at least a subset of the        meteorological data.    -   24. The apparatus of one of 14-23, wherein the meteorological        watch module is adapted to use radiosonde sounding data as at        least a subset of the meteorological data.    -   25. The apparatus of one of 14-24, wherein elapsed time between        obtaining the GNSS data and determining values over the region        for total electron content is no more than about five seconds.    -   26. The apparatus of one of 14-24, wherein elapsed time between        obtaining the GNSS data and determining values over the region        for total electron content and at least one of integrated        precipitable water vapor and tropospheric slant wet delay is no        more than about five seconds.    -   27. A computer program product comprising: a computer usable        medium having computer readable instructions physically embodied        therein, the computer readable instructions when executed by a        processor enabling the processor to perform the method of one of        1-13.    -   28. A computer program comprising a set of instructions which        when loaded in and executed by a processor enable the processor        to perform the method of one of 1-13.        Part 9B: [GNSS Atmospheric Estimation with Ambiguity Fixing]    -   1. A method of processing GNSS signal data to estimate        environmental parameter values, comprising        -   obtaining GNSS data collected at a plurality of stations            distributed over a region from signals received from GNSS            satellites over multiple epochs,        -   obtaining a satellite code bias per satellite,        -   estimating from the GNSS data and the satellite code biases            a set of float ambiguities per station per satellite and            values for atmospheric parameters comprising at least one            of: (i) a total electron content per station per            satellite, (ii) a zenith total delay per station, (iii)a            zenith total delay per station and a set of tropospheric            gradients per station, and (iv) a slant total delay per            station per satellite,        -   fixing the ambiguities, and        -   estimating from the GNSS data and the fixed ambiguities            corrected values for the atmospheric parameters.    -   2. The method of 1, further comprising:        -   obtaining meteorological data for locations within the            region, and        -   determining values over the region for at least one of            integrated precipitable water vapor and tropospheric slant            wet delay, from the corrected values for the atmospheric            parameters.            3. The method of 1 or 2, wherein estimating the set of float            ambiguities and values for atmospheric parameters comprises            applying at least one iterative filter to the GNSS data and            the satellite code biases.    -   4. The method of one of 1-3, wherein estimating the set of float        ambiguities and values for atmospheric parameters comprises        applying a set of factorized filters to the GNSS data and the        satellite code biases.    -   5. The method of one of 1-4, wherein estimating the set of float        ambiguities and values for atmospheric parameters comprises        applying a federated ionospheric filter to the GNSS data and the        satellite code biases.    -   6. The method of one of 1-5, wherein the value of total electron        content per station per satellite comprises a value mapped to        vertical.    -   7. The method of one of 1-6, wherein obtaining meteorological        data comprises obtaining surface meteorological data for        locations within the region.    -   8. The method of one of 1-7, wherein obtaining meteorological        data comprises obtaining radiosonde temperature data for        locations within the region.    -   9. The method of one of 1-8, wherein elapsed time between        obtaining the GNSS data and determining the corrected values for        the atmospheric parameters is no more than about 5 seconds.    -   10. The method of one of 1-9, wherein elapsed time between        obtaining the GNSS data and determining values over the region        for at least one of integrated precipitable water vapor and        tropospheric slant wet delay is no more than about 5 seconds.    -   11. Apparatus for processing GNSS signal data to estimate        environmental parameter values, comprising        -   at least one recursive filter to estimate from satellite            differential code biases and from the GNSS signal data a set            of float ambiguities per station per satellite, and values            for atmospheric parameters comprising at least one of : (i)            a total electron content per station per satellite, (ii) a            zenith total delay per station, (iii) a zenith total delay            per station and a set of tropospheric gradients per            station, (iv) a slant tropospheric total delay per station            per satellite, and (v) a total electron content per station            per satellite.        -   an ambiguity-fixing module to prepare fixed ambiguities from            the float ambiguities, and        -   a corrector module to prepare from the GNSS data and the            fixed ambiguities corrected values for the atmospheric            parameters.    -   12. The apparatus of 11, further comprising:        -   a meteorological watch module to determine values for at            least one of integrated precipitable water vapor and            tropospheric slant wet delay from the corrected values for            the atmospheric parameters and from meteorological data.    -   13. The apparatus of one of 11-12, wherein the at least one        recursive filter comprises a Kalman filter.    -   14. The apparatus of one of 11-13, wherein the at least one        recursive filter comprises a set of factorized filters.    -   15. The apparatus of one of 11-14, wherein the at least one        recursive filter comprises a federated ionospheric filter.    -   16. The apparatus of one of 11-15, wherein the value of total        electron content per station per satellite comprises a value        mapped to vertical.    -   17. The apparatus of one of 11-16, wherein the meteorological        watch module is operative to process meteorological data        comprising surface meteorological data for locations within the        region.    -   18. The apparatus of one of 11-17, wherein the meteorological        watch module is operative to process meteorological data        comprising radiosonde sounding data for locations within the        region.

19. The apparatus of one of 11-18, wherein elapsed time betweenobtaining the GNSS data and determining the corrected values over theregion for the atmospheric parameters is no more than about 5 seconds.

-   -   20. The apparatus of one of 11-19, wherein elapsed time between        obtaining the GNSS data and determining the corrected values        over the region for at least one of integrated precipitable        water vapor and tropospheric slant wet delay is no more than        about five seconds    -   21. A computer program product comprising: a computer usable        medium having computer readable instructions physically embodied        therein, the computer readable instructions when executed by a        processor enabling the processor to perform the method of one of        1-10.    -   22. A computer program comprising a set of instructions which        when loaded in and executed by a processor enable the processor        to perform the method of one of 1-10.

1. A method of processing GNSS signal data to estimate environmentalparameter values, comprising obtaining GNSS data collected at aplurality of stations distributed over a region from signals receivedfrom GNSS satellites over multiple epochs, obtaining a satellite codebias per satellite, estimating from the GNSS data and the satellite codebiases a set of float ambiguities per station per satellite and valuesfor atmospheric parameters comprising at least one of: (i) a totalelectron content per station per satellite, (ii) a zenith total delayper station, (iii) a zenith total delay per station and a set oftropospheric gradients per station, and (iv) a slant total delay perstation per satellite, fixing the ambiguities, and estimating from theGNSS data and the fixed ambiguities corrected values for the atmosphericparameters.
 2. The method of claim 1, further comprising: obtainingmeteorological data for locations within the region, and determiningvalues over the region for at least one of integrated precipitable watervapor and tropospheric slant wet delay, from the corrected values forthe atmospheric parameters and from the meteorological data.
 3. Themethod of claim 1, wherein estimating the set of float ambiguities andvalues for atmospheric parameters comprises applying at least oneiterative filter to the GNSS data and the satellite code biases.
 4. Themethod of claim 1, wherein estimating the set of float ambiguities andvalues for atmospheric parameters comprises applying a set of factorizedfilters to the GNSS data and the satellite code biases.
 5. The method ofclaim 1, wherein estimating the set of float ambiguities and values foratmospheric parameters comprises applying a federated ionospheric filterto the GNSS data and the satellite code biases.
 6. The method of claim1, wherein the value of total electron content per station per satellitecomprises a value mapped to vertical.
 7. The method of claim 1, whereinobtaining meteorological data comprises obtaining surface meteorologicaldata for locations within the region.
 8. The method of claim 1, whereinobtaining meteorological data comprises obtaining radiosonde temperaturedata for locations within the region.
 9. The method of claim 1, whereinelapsed time between obtaining the GNSS data and determining thecorrected values for the atmospheric parameters is no more than about 5seconds.
 10. The method of claim 1, wherein elapsed time betweenobtaining the GNSS data and determining values over the region for atleast one of integrated precipitable water vapor and tropospheric slantwet delay is no more than about 5 seconds.
 11. Apparatus for processingGNSS signal data to estimate environmental parameter values, comprisingat least one recursive filter to estimate from satellite differentialcode biases and from the GNSS signal data a set of float ambiguities perstation per satellite, and values for atmospheric parameters comprisingat least one of : (i) a total electron content per station persatellite, (ii) a zenith total delay per station, (iii) a zenith totaldelay per station and a set of tropospheric gradients per station, (iv)a slant tropospheric total delay per station per satellite, and (v) atotal electron content per station per satellite. an ambiguity-fixingmodule to prepare fixed ambiguities from the float ambiguities, and acorrector module to prepare from the GNSS data and the fixed ambiguitiescorrected values for the atmospheric parameters.
 12. The apparatus ofclaim 11, further comprising: a meteorological watch module to determinevalues for at least one of integrated precipitable water vapor andtropospheric slant wet delay from the corrected values for theatmospheric parameters and from meteorological data.
 13. The apparatusof claim 11, wherein the at least one recursive filter comprises aKalman filter.
 14. The apparatus of claim 11, wherein the at least onerecursive filter comprises a set of factorized filters.
 15. Theapparatus of claim 11, wherein the at least one recursive filtercomprises a federated ionospheric filter.
 16. The apparatus of claim 11,wherein the value of total electron content per station per satellitecomprises a value mapped to vertical.
 17. The apparatus of claim 11,wherein the meteorological watch module is operative to processmeteorological data comprising surface meteorological data for locationswithin the region.
 18. The apparatus of claim 11, wherein themeteorological watch module is operative to process meteorological datacomprising radiosonde sounding data for locations within the region. 19.The apparatus of claim 11, wherein elapsed time between obtaining theGNSS data and determining the corrected values over the region for theatmospheric parameters is no more than about 5 seconds.
 20. Theapparatus of claim 11, wherein elapsed time between obtaining the GNSSdata and determining values over the region for integrated precipitablewater vapor and ionospheric slant wet delay is no more than about 5seconds.
 21. A computer program product comprising: a computer usablemedium having computer readable instructions physically embodiedtherein, the computer readable instructions when executed by a processorenabling the processor to perform the method of claim 1.